Setup chunk
Setup reticulate
knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")
[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/pallium_evo"
Load libraries
library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)
[1] FALSE
use_python(Sys.which("python"))
py_config()
python: /home/tpires/bin/miniconda3/bin/python
libpython: /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome: /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version: 3.8.3 (default, May 19 2020, 18:47:26) [GCC 7.3.0]
numpy: /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version: 1.18.5
NOTE: Python version was forced by RETICULATE_PYTHON
Set colours for cell types and regions
library(Seurat)
Registered S3 method overwritten by 'spatstat.geom':
method from
print.boxx cli
Attaching SeuratObject
library(patchwork)
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(tidyverse)
── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5 ✓ purrr 0.3.4
✓ tibble 3.1.6 ✓ stringr 1.4.0
✓ tidyr 1.1.4 ✓ forcats 0.5.1
✓ readr 2.1.0
── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(ggplot2)
library(parallel)
library(doParallel)
Loading required package: foreach
Attaching package: ‘foreach’
The following objects are masked from ‘package:purrr’:
accumulate, when
Loading required package: iterators
library(foreach)
library(ggdendro)
library(presto)
Loading required package: Rcpp
Loading required package: data.table
data.table 1.14.2 using 72 threads (see ?getDTthreads). Latest news: r-datatable.com
Attaching package: ‘data.table’
The following object is masked from ‘package:purrr’:
transpose
The following objects are masked from ‘package:dplyr’:
between, first, last
library(scatterpie)
Load Div-seq
meta = read.csv("data/annotations/axolotl_all_umeta.csv",
header = T, row.names = 1)
cols_cc = c(
#epen
"#12400c", "#2d6624","#1d4f15", "#174711", "#2d6624", "#3d7f33", "#3b7b30", "#468b3b", "#4f9843","#5dae50", "#66bb58", "#72cd64", "#306a26", "#78d669", "#81e472",
#gaba
"#700209", "#75090e","#7a0f13", "#801517", "#851a1b", "#8a1f1f", "#902423", "#952927", "#9a2d2c","#a03230", "#a53634", "#aa3a39", "#b03f3d","#b54342", "#ba4846", "#c04c4b", "#c5504f", "#ca5554", "#d05959", "#d55e5e","#73050c", "#780c11","#8d2221", "#982b2a","#a23432", "#a83837", "#b2413f", "#b84544", "#bd4a49", "#c85352", #"#cd5756",
#glut
"#054674", "#134d7b","#1d5481", "#265a88", "#2e618e", "#73a4cb", "#366995", "#3e709c", "#4677a2","#4d7ea9", "#5586b0", "#5c8db7", "#6495bd","#6b9cc4", "#7bacd2", "#8ebfe4", "#96c7eb", "#9ecff2", "#18507e", "#18507e","#2a5e8b", "#497ba6","#5889b3", "#6fa0c8","#7fafd6", "#6091ba", "#5182ac", "#3a6c98", "#a6d7f9",
#npc
"#ffb120", "#feb72a","#fdbc34", "#fcc13d", "#fbc745", "#facc4e", "#f9d156", "#f8d65f", "#f8da68","#f7df70", "#f7e479", "#f7e882", "#f7ed8a", "#f7f193", "#eca319"
)
ccnames = unique(sort(meta$cellclusters))
names(cols_cc) = c(ccnames[grepl("epen", ccnames)], ccnames[grepl("GABA", ccnames)],ccnames[grepl("glut", ccnames)],ccnames[grepl("npc", ccnames)])
cols_cc = c(cols_cc, "microglia_8" = "#E6530D",
"oligodendrocyte_15" = "#E43D88", "oligodendrocyte_10" = "#F662A5",
"endothelial_11" = "#712166", "endothelial_12" = "#B0279D",
"endothelial_14" = "#BE5AB0")
reg_cols = c("other/unknown_pred" = "#C7CCC7",
"medial" = "#52168D", "medial_pred" = "#661CB0",
"dorsal" = "#C56007", "dorsal_pred" = "#ED7307",
"lateral" = "#118392", "lateral_pred" = "#16A3B6")
reg_cols_simp = c("medial" = "#52168D", "dorsal" = "#C56007", "lateral" = "#118392")
Save necessary data
```r
meta = read.csv(\data/annotations/axolotl_all_umeta.csv\,
header = T, row.names = 1)
cols_cc = c(
#epen
\#12400c\, \#2d6624\,\#1d4f15\, \#174711\, \#2d6624\, \#3d7f33\, \#3b7b30\, \#468b3b\, \#4f9843\,\#5dae50\, \#66bb58\, \#72cd64\, \#306a26\, \#78d669\, \#81e472\,
#gaba
\#700209\, \#75090e\,\#7a0f13\, \#801517\, \#851a1b\, \#8a1f1f\, \#902423\, \#952927\, \#9a2d2c\,\#a03230\, \#a53634\, \#aa3a39\, \#b03f3d\,\#b54342\, \#ba4846\, \#c04c4b\, \#c5504f\, \#ca5554\, \#d05959\, \#d55e5e\,\#73050c\, \#780c11\,\#8d2221\, \#982b2a\,\#a23432\, \#a83837\, \#b2413f\, \#b84544\, \#bd4a49\, \#c85352\, #\#cd5756\,
#glut
\#054674\, \#134d7b\,\#1d5481\, \#265a88\, \#2e618e\, \#73a4cb\, \#366995\, \#3e709c\, \#4677a2\,\#4d7ea9\, \#5586b0\, \#5c8db7\, \#6495bd\,\#6b9cc4\, \#7bacd2\, \#8ebfe4\, \#96c7eb\, \#9ecff2\, \#18507e\, \#18507e\,\#2a5e8b\, \#497ba6\,\#5889b3\, \#6fa0c8\,\#7fafd6\, \#6091ba\, \#5182ac\, \#3a6c98\, \#a6d7f9\,
#npc
\#ffb120\, \#feb72a\,\#fdbc34\, \#fcc13d\, \#fbc745\, \#facc4e\, \#f9d156\, \#f8d65f\, \#f8da68\,\#f7df70\, \#f7e479\, \#f7e882\, \#f7ed8a\, \#f7f193\, \#eca319\
)
ccnames = unique(sort(meta$cellclusters))
names(cols_cc) = c(ccnames[grepl(\epen\, ccnames)], ccnames[grepl(\GABA\, ccnames)],ccnames[grepl(\glut\, ccnames)],ccnames[grepl(\npc\, ccnames)])
cols_cc = c(cols_cc, \microglia_8\ = \#E6530D\,
\oligodendrocyte_15\ = \#E43D88\, \oligodendrocyte_10\ = \#F662A5\,
\endothelial_11\ = \#712166\, \endothelial_12\ = \#B0279D\,
\endothelial_14\ = \#BE5AB0\)
reg_cols = c(\other/unknown_pred\ = \#C7CCC7\,
\medial\ = \#52168D\, \medial_pred\ = \#661CB0\,
\dorsal\ = \#C56007\, \dorsal_pred\ = \#ED7307\,
\lateral\ = \#118392\, \lateral_pred\ = \#16A3B6\)
reg_cols_simp = c(\medial\ = \#52168D\, \dorsal\ = \#C56007\, \lateral\ = \#118392\)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
# Region predictions
Load results from models
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxub3V0bmFtZSA9IFxcZGF0YS9wcm9jZXNzZWQvYXhvbG90bF9wYXJ0cy9kaXZfcmVnaW9uc19kYXRhLm10eFxcXG5vdXRuYW1lMiA9IFxcZGF0YS9wcm9jZXNzZWQvYXhvbG90bF9wYXJ0cy9kaXZfcmVnaW9uc19jb3VudHMubXR4XFxcbmlmKCFmaWxlLmV4aXN0cyhvdXRuYW1lKSl7XG4gIG1ldGFkYXRhX2F4ID0gZGl2X2RhdEBtZXRhLmRhdGFbLGMoXFxzYW1wbGVcXCwgXFxiYXRjaFxcLCBcXGhpZ2hsZXZlbF9jbHVzdGVyc1xcKV1cbiAgd3JpdGUuY3N2KG1ldGFkYXRhX2F4LCBmaWxlID0gXFxkYXRhL3Byb2Nlc3NlZC9heG9sb3RsX3BhcnRzL2Rpdl9yZWdpb25zX21ldGEuY3N2XFwsIFxuICAgICAgICAgICAgY29sLm5hbWVzID0gVCwgcm93Lm5hbWVzID0gVCwgcXVvdGUgPSBGKVxuICB3cml0ZS5jc3YoZGl2X2RhdEBhc3NheXMkUk5BQG1ldGEuZmVhdHVyZXMsIFxuICAgICAgICAgICAgZmlsZSA9IFxcZGF0YS9wcm9jZXNzZWQvYXhvbG90bF9wYXJ0cy9kaXZfcmVnaW9uc19nZW5lcy5jc3ZcXCwgXG4gICAgICAgICAgICBjb2wubmFtZXMgPSBULCByb3cubmFtZXMgPSBULCBxdW90ZSA9IEYpXG4gIFxuICBNYXRyaXg6OndyaXRlTU0oTWF0cml4Ojp0KGRpdl9kYXRAYXNzYXlzJFJOQUBkYXRhKSwgb3V0bmFtZSlcbiAgTWF0cml4Ojp3cml0ZU1NKE1hdHJpeDo6dChkaXZfZGF0QGFzc2F5cyRSTkFAY291bnRzKSwgb3V0bmFtZTIpXG59XG5gYGBcbmBgYCJ9 -->
```r
```r
outname = \data/processed/axolotl_parts/div_regions_data.mtx\
outname2 = \data/processed/axolotl_parts/div_regions_counts.mtx\
if(!file.exists(outname)){
metadata_ax = div_dat@meta.data[,c(\sample\, \batch\, \highlevel_clusters\)]
write.csv(metadata_ax, file = \data/processed/axolotl_parts/div_regions_meta.csv\,
col.names = T, row.names = T, quote = F)
write.csv(div_dat@assays$RNA@meta.features,
file = \data/processed/axolotl_parts/div_regions_genes.csv\,
col.names = T, row.names = T, quote = F)
Matrix::writeMM(Matrix::t(div_dat@assays$RNA@data), outname)
Matrix::writeMM(Matrix::t(div_dat@assays$RNA@counts), outname2)
}
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Confusion matrices
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYWxsX3ByZWRfZiA9IGxpc3QuZmlsZXMoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHBhdHRlcm4gPSBcXHByZWRzX1xcKVxuYWxsX3ByZWRfZiA9IGFsbF9wcmVkX2ZbZ3JlcGwoXFxsclxcLCBhbGxfcHJlZF9mKSAmICFncmVwbChcXENUXFwsIGFsbF9wcmVkX2YpXVxuIyB0ZXN0IHNldCByZXN1bHRzXG50ZXN0X3JlcyA9IGxhcHBseShhbGxfcHJlZF9mW2dyZXBsKFxcYXguY3N2XFwsIGFsbF9wcmVkX2YpXSwgXG4gICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKHRlc3RfcmVzKSA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQoYWxsX3ByZWRfZltncmVwbChcXGF4LmNzdlxcLCBhbGxfcHJlZF9mKV0sIFxcLlxcLCBmaXhlZCA9IFQpLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSB4WzFdKSlcblxuIyBEaXYtc2VxIHJlc3VsdHNcbmRpdl9yZXMgPSBsYXBwbHkoYWxsX3ByZWRfZltncmVwbChcXGRpdlxcLCBhbGxfcHJlZF9mKV0sIFxuICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKGRpdl9yZXMpID0gdW5saXN0KGxhcHBseShzdHJzcGxpdChhbGxfcHJlZF9mW2dyZXBsKFxcZGl2XFwsIGFsbF9wcmVkX2YpXSwgXFwuXFwsIGZpeGVkID0gVCksXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgeFsxXSkpXG5gYGBcbmBgYCJ9 -->
```r
```r
all_pred_f = list.files(\results/Div-seq/\, pattern = \preds_\)
all_pred_f = all_pred_f[grepl(\lr\, all_pred_f) & !grepl(\CT\, all_pred_f)]
# test set results
test_res = lapply(all_pred_f[grepl(\ax.csv\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(test_res) = unlist(lapply(strsplit(all_pred_f[grepl(\ax.csv\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
# Div-seq results
div_res = lapply(all_pred_f[grepl(\div\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(div_res) = unlist(lapply(strsplit(all_pred_f[grepl(\div\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Compare all with neu/dif (test sets are different so overlap is not complete)
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYWxsX3ByZWRfZiA9IGxpc3QuZmlsZXMoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHBhdHRlcm4gPSBcXHByZWRzX1xcKVxuYWxsX3ByZWRfZiA9IGFsbF9wcmVkX2ZbZ3JlcGwoXFxsclxcLCBhbGxfcHJlZF9mKSAmICFncmVwbChcXENUXFwsIGFsbF9wcmVkX2YpXVxuIyB0ZXN0IHNldCByZXN1bHRzXG50ZXN0X3JlcyA9IGxhcHBseShhbGxfcHJlZF9mW2dyZXBsKFxcYXguY3N2XFwsIGFsbF9wcmVkX2YpXSwgXG4gICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKHRlc3RfcmVzKSA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQoYWxsX3ByZWRfZltncmVwbChcXGF4LmNzdlxcLCBhbGxfcHJlZF9mKV0sIFxcLlxcLCBmaXhlZCA9IFQpLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSB4WzFdKSlcblxuIyBEaXYtc2VxIHJlc3VsdHNcbmRpdl9yZXMgPSBsYXBwbHkoYWxsX3ByZWRfZltncmVwbChcXGRpdlxcLCBhbGxfcHJlZF9mKV0sIFxuICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKGRpdl9yZXMpID0gdW5saXN0KGxhcHBseShzdHJzcGxpdChhbGxfcHJlZF9mW2dyZXBsKFxcZGl2XFwsIGFsbF9wcmVkX2YpXSwgXFwuXFwsIGZpeGVkID0gVCksXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgeFsxXSkpXG5gYGBcbmBgYCJ9 -->
```r
```r
all_pred_f = list.files(\results/Div-seq/\, pattern = \preds_\)
all_pred_f = all_pred_f[grepl(\lr\, all_pred_f) & !grepl(\CT\, all_pred_f)]
# test set results
test_res = lapply(all_pred_f[grepl(\ax.csv\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(test_res) = unlist(lapply(strsplit(all_pred_f[grepl(\ax.csv\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
# Div-seq results
div_res = lapply(all_pred_f[grepl(\div\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(div_res) = unlist(lapply(strsplit(all_pred_f[grepl(\div\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Plot UMAP with regions
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxudGVzdF9sYWJzX2wgPSBsaXN0KClcbmZvcihnIGluIGMoXFxhbGxcXCwgXFxkaWZcXCwgXFxuZXVcXCkpe1xuICBuMSA9IG5hbWVzKHRlc3RfcmVzKVtncmVwbChnLCBuYW1lcyh0ZXN0X3JlcykpICYgZ3JlcGwoXFxsclxcLCBuYW1lcyh0ZXN0X3JlcykpXVxuICB0ZXN0X2xhYnNfbFtbZ11dID0gdGVzdF9yZXNbW24xXV1cbiAgcHJpbnQoZylcbiAgcHJpbnQodGFibGUodGVzdF9sYWJzX2xbW2ddXSR5X3Rlc3QsIHRlc3RfbGFic19sW1tnXV0kcHJlZF9scikpXG59XG5gYGBcbmBgYCJ9 -->
```r
```r
test_labs_l = list()
for(g in c(\all\, \dif\, \neu\)){
n1 = names(test_res)[grepl(g, names(test_res)) & grepl(\lr\, names(test_res))]
test_labs_l[[g]] = test_res[[n1]]
print(g)
print(table(test_labs_l[[g]]$y_test, test_labs_l[[g]]$pred_lr))
}
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiWzFdIFxcYWxsXFxcbiAgICAgICAgIFxuICAgICAgICAgIGRvcnNhbCBsYXRlcmFsIG1lZGlhbFxuICBkb3JzYWwgICAgMjAxNyAgICAgIDc3ICAgICA4MVxuICBsYXRlcmFsICAgICA2OSAgICAzMDIwICAgICAyMlxuICBtZWRpYWwgICAgICA2MyAgICAgIDY1ICAgMzc1M1xuWzFdIFxcZGlmXFxcbiAgICAgICAgIFxuICAgICAgICAgIGRvcnNhbCBsYXRlcmFsIG1lZGlhbFxuICBkb3JzYWwgICAgMTY4MyAgICAgIDU5ICAgICA5NFxuICBsYXRlcmFsICAgICA2MSAgICAyNjU0ICAgICAyMVxuICBtZWRpYWwgICAgICA1NyAgICAgIDU0ICAgMzQzNVxuWzFdIFxcbmV1XFxcbiAgICAgICAgIFxuICAgICAgICAgIGRvcnNhbCBsYXRlcmFsIG1lZGlhbFxuICBkb3JzYWwgICAgMTczNiAgICAgIDY1ICAgICA5MVxuICBsYXRlcmFsICAgICA3NyAgICAyNzYzICAgICAzNFxuICBtZWRpYWwgICAgICA1MCAgICAgIDUxICAgMzUwOVxuIn0= -->
[1]
dorsal lateral medial
dorsal 2017 77 81 lateral 69 3020 22 medial 63 65 3753 [1]
dorsal lateral medial
dorsal 1683 59 94 lateral 61 2654 21 medial 57 54 3435 [1]
dorsal lateral medial
dorsal 1736 65 91 lateral 77 2763 34 medial 50 51 3509
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGl2X2xhYnNfbCA9IGxpc3QoKVxuZm9yKGcgaW4gYyhcXGFsbFxcLCBcXGRpZlxcLCBcXG5ldVxcKSl7XG4gIG4xID0gbmFtZXMoZGl2X3JlcylbZ3JlcGwoZywgbmFtZXMoZGl2X3JlcykpICYgZ3JlcGwoXFxsclxcLCBuYW1lcyhkaXZfcmVzKSldXG4gIGRpdl9sYWJzX2xbW2ddXSA9IGRpdl9yZXNbW24xXV1cbn0gXG5gYGBcbmBgYCJ9 -->
```r
```r
div_labs_l = list()
for(g in c(\all\, \dif\, \neu\)){
n1 = names(div_res)[grepl(g, names(div_res)) & grepl(\lr\, names(div_res))]
div_labs_l[[g]] = div_res[[n1]]
}
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Plot changes in proportion
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxudGVzdF9hbGxfbmV1ID0gbWVyZ2UodGVzdF9sYWJzX2wkYWxsLCB0ZXN0X2xhYnNfbCRuZXUsIGJ5ID0gMSlcbnRhYmxlKHRlc3RfYWxsX25ldSRwcmVkX2xyLngsIHRlc3RfYWxsX25ldSRwcmVkX2xyLnkpXG5gYGBcbmBgYCJ9 -->
```r
```r
test_all_neu = merge(test_labs_l$all, test_labs_l$neu, by = 1)
table(test_all_neu$pred_lr.x, test_all_neu$pred_lr.y)
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgXG4gICAgICAgICAgZG9yc2FsIGxhdGVyYWwgbWVkaWFsXG4gIGRvcnNhbCAgICAgNDY2ICAgICAgIDMgICAgIDEzXG4gIGxhdGVyYWwgICAgICA4ICAgICA5NjEgICAgIDEwXG4gIG1lZGlhbCAgICAgICA1ICAgICAgIDQgICAxMzM4XG4ifQ== -->
dorsal lateral medial
dorsal 466 3 13 lateral 8 961 10 medial 5 4 1338
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxudGVzdF9hbGxfZXAgPSBtZXJnZSh0ZXN0X2xhYnNfbCRhbGwsIHRlc3RfbGFic19sJGRpZiwgYnkgPSAxKVxudGFibGUodGVzdF9hbGxfZXAkcHJlZF9sci54LCB0ZXN0X2FsbF9lcCRwcmVkX2xyLnkpXG5gYGBcbmBgYCJ9 -->
```r
```r
test_all_ep = merge(test_labs_l$all, test_labs_l$dif, by = 1)
table(test_all_ep$pred_lr.x, test_all_ep$pred_lr.y)
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgXG4gICAgICAgICAgZG9yc2FsIGxhdGVyYWwgbWVkaWFsXG4gIGRvcnNhbCAgICAgNDczICAgICAgIDUgICAgIDExXG4gIGxhdGVyYWwgICAgICAzICAgICA5MzIgICAgIDExXG4gIG1lZGlhbCAgICAgICA2ICAgICAgIDYgICAxMzM0XG4ifQ== -->
dorsal lateral medial
dorsal 473 5 11 lateral 3 932 11 medial 6 6 1334
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGl2X2FsbF9uZXUgPSBtZXJnZShkaXZfbGFic19sJGFsbCwgZGl2X2xhYnNfbCRuZXUsIGJ5ID0gMSlcbnRhYmxlKGRpdl9hbGxfbmV1JHByZWRfbHIueCwgZGl2X2FsbF9uZXUkcHJlZF9sci55KVxuYGBgXG5gYGAifQ== -->
```r
```r
div_all_neu = merge(div_labs_l$all, div_labs_l$neu, by = 1)
table(div_all_neu$pred_lr.x, div_all_neu$pred_lr.y)
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgXG4gICAgICAgICAgZG9yc2FsIGxhdGVyYWwgbWVkaWFsXG4gIGRvcnNhbCAgICA4NjgxICAgICA0NjAgICAgMzI2XG4gIGxhdGVyYWwgICAgNDQzICAgIDY2MDIgICAgIDk2XG4gIG1lZGlhbCAgICAgMTg0ICAgICAxMTEgICAyMzk3XG4ifQ== -->
dorsal lateral medial
dorsal 8681 460 326 lateral 443 6602 96 medial 184 111 2397
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGl2X2FsbF9lcCA9IG1lcmdlKGRpdl9sYWJzX2wkYWxsLCBkaXZfbGFic19sJGRpZiwgYnkgPSAxKVxudGFibGUoZGl2X2FsbF9lcCRwcmVkX2xyLngsIGRpdl9hbGxfZXAkcHJlZF9sci55KVxuYGBgXG5gYGAifQ== -->
```r
```r
div_all_ep = merge(div_labs_l$all, div_labs_l$dif, by = 1)
table(div_all_ep$pred_lr.x, div_all_ep$pred_lr.y)
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgXG4gICAgICAgICAgZG9yc2FsIGxhdGVyYWwgbWVkaWFsXG4gIGRvcnNhbCAgICA4NzQ2ICAgICAzNjQgICAgMzU3XG4gIGxhdGVyYWwgICAgNzEwICAgIDYyOTUgICAgMTM2XG4gIG1lZGlhbCAgICAgMjIxICAgICAxMDYgICAyMzY1XG4ifQ== -->
dorsal lateral medial
dorsal 8746 364 357 lateral 710 6295 136 medial 221 106 2365
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
# Cell Type predictions
Load results from models
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGl2X2RhdCRwcmVkX3JlZ3MgPSBkaXZfcmVzJHByZWRzX2xyX2Rpdl9hbGwkcHJlZF9sclxuRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxcaGlnaGxldmVsX2NsdXN0ZXJzXFwsIGxhYmVsID0gVClcbmBgYFxuYGBgIn0= -->
```r
```r
div_dat$pred_regs = div_res$preds_lr_div_all$pred_lr
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \highlevel_clusters\, label = T)
<!-- rnb-source-end -->
<!-- rnb-plot-begin -->
<img src="data:image/png;base64," />
<!-- rnb-plot-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxccHJlZF9yZWdzXFwpXG5gYGBcbmBgYCJ9 -->
```r
```r
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \pred_regs\)
<!-- rnb-source-end -->
<!-- rnb-plot-begin -->
<img src="data:image/png;base64," />
<!-- rnb-plot-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuRmVhdHVyZVBsb3QoZGl2X2RhdCwgcmVkdWN0aW9uID0gXFx1bWFwX2hhcm1vbnlcXCwgZmVhdHVyZXMgPSBjKFxcU0ZSUDFcXCwgXFxXTlQ4QlxcLCBcXFVOQzVDXFwpKVxuXG5yZWdfY29sc19zaW1wID0gYyhcXG1lZGlhbFxcID0gXFwjNTIxNjhEXFwsIFxcZG9yc2FsXFwgPSBcXCNDNTYwMDdcXCwgXFxsYXRlcmFsXFwgPSBcXCMxMTgzOTJcXClcblxucGx0ID0gRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxccHJlZF9yZWdzXFwsIHNodWZmbGUgPSBUKStcbiAgc2NhbGVfY29sb3VyX21hbnVhbCh2YWx1ZXMgPSByZWdfY29sc19zaW1wKStcbiAgdGhlbWUoYXNwZWN0LnJhdGlvID0gMSxcbiAgICAgICAgYXhpcy5saW5lID0gZWxlbWVudF9ibGFuaygpLFxuICAgICAgICBheGlzLnRpY2tzID0gZWxlbWVudF9ibGFuaygpLFxuICAgICAgICBheGlzLnRleHQgPSBlbGVtZW50X2JsYW5rKCksXG4gICAgICAgIGF4aXMudGl0bGUgPSBlbGVtZW50X2JsYW5rKCksXG4gICAgICAgIHBsb3QudGl0bGUgPSBlbGVtZW50X2JsYW5rKCksXG4gICAgICAgIGxlZ2VuZC5wb3NpdGlvbiA9IFxcbm9uZVxcKVxuXG5wZGYoXFxyZXN1bHRzL0Rpdi1zZXEvRGl2c2VxX3Jlc3VsdHMvZGl2c2VxX3VtYXBfcmVnaW9ucy5wZGZcXCwgaGVpZ2h0ID0gMi43LCB3aWR0aCA9IDIuNylcbmBgYFxuYGBgclxucHJpbnQocGx0KVxuZGV2Lm9mZigpXG5gYGBcbmBgYCJ9 -->
```r
```r
FeaturePlot(div_dat, reduction = \umap_harmony\, features = c(\SFRP1\, \WNT8B\, \UNC5C\))
reg_cols_simp = c(\medial\ = \#52168D\, \dorsal\ = \#C56007\, \lateral\ = \#118392\)
plt = DimPlot(div_dat, reduction = \umap_harmony\, group.by = \pred_regs\, shuffle = T)+
scale_colour_manual(values = reg_cols_simp)+
theme(aspect.ratio = 1,
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
plot.title = element_blank(),
legend.position = \none\)
pdf(\results/Div-seq/Divseq_results/divseq_umap_regions.pdf\, height = 2.7, width = 2.7)
print(plt)
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoicG5nIFxuICAyIFxuIn0= -->
png 2
<!-- rnb-output-end -->
<!-- rnb-plot-begin -->
<img src="data:image/png;base64," />
<!-- rnb-plot-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Plot predictions
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuXG50YWJfY3RfcmVnID0gdGFibGUocGFzdGUwKGRpdl9kYXQkc2FtcGxlLCBcXC4uXFwsIGRpdl9kYXQkaGlnaF9sZXZlbF9hbm5vKSwgZGl2X2RhdCRwcmVkX3JlZ3MpXG50YWJfY3RfcmVnID0gZGF0YS5mcmFtZSh0YWJfY3RfcmVnKVxuXG50YWJfY3RfcmVnJHRpbWUgPSB1bmxpc3QobGFwcGx5KHN0cnNwbGl0KGFzLmNoYXJhY3Rlcih0YWJfY3RfcmVnJFZhcjEpLCBcXF9cXCksIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSBwYXN0ZTAoeFsxXSwgXFxfd3BpXFwpKSlcbnRhYl9jdF9yZWckY3QgPSB1bmxpc3QobGFwcGx5KHN0cnNwbGl0KGFzLmNoYXJhY3Rlcih0YWJfY3RfcmVnJFZhcjEpLCBcXC4uXFwsIGZpeGVkID0gVCksIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSB4WzJdKSlcblxuZm9yKHIgaW4gdW5pcXVlKHRhYl9jdF9yZWckdGltZSkpe1xuICBjYyA9IHRhYl9jdF9yZWckdGltZT09clxuICB0YWJfY3RfcmVnJEZyZXFbY2NdID0gdGFiX2N0X3JlZyRGcmVxW2NjXS9zdW0odGFiX2N0X3JlZyRGcmVxW2NjXSlcbn1cblxuZm9yKHIgaW4gdW5pcXVlKHRhYl9jdF9yZWckY3QpKXtcbiAgY2MgPSB0YWJfY3RfcmVnJGN0PT1yXG4gIHRhYl9jdF9yZWckRnJlcVtjY10gPSBzY2FsZXM6OnJlc2NhbGUodGFiX2N0X3JlZyRGcmVxW2NjXSwgYygwLDEpKVxufVxuXG50YWJfY3RfcmVnJGN0ID0gZmFjdG9yKHRhYl9jdF9yZWckY3QsIGxldmVscyA9IGMoXFxtaWNyb2dsaWFcXCwgXFxlbmRvXFwsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFxcZXBlbmR5bWFsX2FcXCwgXFxlcGVuZHltYWxfcVxcLCBcXE5QQ1xcLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFxcbmV1cm9uYWxcXCkpXG50YWJfY3RfcmVnJHRpbWUgPSBmYWN0b3IodGFiX2N0X3JlZyR0aW1lLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICBsZXZlbHMgPSByZXYoYyhcXDFfd3BpXFwsIFxcMl93cGlcXCwgXFw0X3dwaVxcLCBcXDZfd3BpXFwsIFxcOF93cGlcXCwgXFwxMl93cGlcXCkpKVxudGFiX2N0X3JlZyRWYXIyID0gZmFjdG9yKHRhYl9jdF9yZWckVmFyMiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgbGV2ZWxzID0gYyhcXG1lZGlhbFxcLFxcZG9yc2FsXFwsXFxsYXRlcmFsXFwpKVxuXG5wbHQgPSBnZ3Bsb3QodGFiX2N0X3JlZywgYWVzKHggPSBWYXIyLCB5ID0gdGltZSwgXG4gICAgICAgICAgICAgICAgICAgICAgIHNpemUgPSBGcmVxLCBhbHBoYSA9IEZyZXEsIGNvbG91ciA9IFZhcjIpKStcbiAgZmFjZXRfd3JhcCh+Y3QsIG5jb2wgPSA2KStcbiAgc2NhbGVfY29sb3VyX21hbnVhbCh2YWx1ZXMgPSByZWdfY29sc19zaW1wKStcbiAgZ2VvbV9wb2ludCgpK1xuICBsYWJzKHggPSBcXFJlZ2lvblxcLCB5ID0gXFxUaW1lXFwsIGNvbG91ciA9IFxcUmVnaW9uXFwsIFxuICAgICAgIHNpemUgPSBcXE5vcm0uIHByb3BvcnRpb25cXCwgYWxwaGEgPSBcXE5vcm0uIHByb3BvcnRpb25cXCkrXG4gIHRoZW1lX2J3KCkrXG4gIHRoZW1lKGF4aXMudGV4dCA9IGVsZW1lbnRfdGV4dChjb2xvdXIgPSBcXGJsYWNrXFwsIHNpemUgPSA3KSxcbiAgICAgICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSAzMCwgaGp1c3QgPSAxLCB2anVzdCA9IDEpLFxuICAgICAgICBheGlzLnRpdGxlID0gZWxlbWVudF90ZXh0KHNpemUgPSA3KSxcbiAgICAgICAgc3RyaXAudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNyksXG4gICAgICAgIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplID0gNyksXG4gICAgICAgIGxlZ2VuZC50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemUgPSA2LjUpKVxuXG5wZGYoXFxyZXN1bHRzL0Rpdi1zZXEvRGl2c2VxX3Jlc3VsdHMvcHJvcG9ydGlvbl9yZWdpb25fbWFqb3JjdF90aW1lLnBkZlxcLCBoZWlnaHQgPSAyLjcsIHdpZHRoID0gNy41KVxucHJpbnQocGx0KVxuZGV2Lm9mZigpXG5gYGBcbmBgYCJ9 -->
```r
```r
tab_ct_reg = table(paste0(div_dat$sample, \..\, div_dat$high_level_anno), div_dat$pred_regs)
tab_ct_reg = data.frame(tab_ct_reg)
tab_ct_reg$time = unlist(lapply(strsplit(as.character(tab_ct_reg$Var1), \_\),
function(x) paste0(x[1], \_wpi\)))
tab_ct_reg$ct = unlist(lapply(strsplit(as.character(tab_ct_reg$Var1), \..\, fixed = T),
function(x) x[2]))
for(r in unique(tab_ct_reg$time)){
cc = tab_ct_reg$time==r
tab_ct_reg$Freq[cc] = tab_ct_reg$Freq[cc]/sum(tab_ct_reg$Freq[cc])
}
for(r in unique(tab_ct_reg$ct)){
cc = tab_ct_reg$ct==r
tab_ct_reg$Freq[cc] = scales::rescale(tab_ct_reg$Freq[cc], c(0,1))
}
tab_ct_reg$ct = factor(tab_ct_reg$ct, levels = c(\microglia\, \endo\,
\ependymal_a\, \ependymal_q\, \NPC\,
\neuronal\))
tab_ct_reg$time = factor(tab_ct_reg$time,
levels = rev(c(\1_wpi\, \2_wpi\, \4_wpi\, \6_wpi\, \8_wpi\, \12_wpi\)))
tab_ct_reg$Var2 = factor(tab_ct_reg$Var2,
levels = c(\medial\,\dorsal\,\lateral\))
plt = ggplot(tab_ct_reg, aes(x = Var2, y = time,
size = Freq, alpha = Freq, colour = Var2))+
facet_wrap(~ct, ncol = 6)+
scale_colour_manual(values = reg_cols_simp)+
geom_point()+
labs(x = \Region\, y = \Time\, colour = \Region\,
size = \Norm. proportion\, alpha = \Norm. proportion\)+
theme_bw()+
theme(axis.text = element_text(colour = \black\, size = 7),
axis.text.x = element_text(angle = 30, hjust = 1, vjust = 1),
axis.title = element_text(size = 7),
strip.text = element_text(size = 7),
legend.title = element_text(size = 7),
legend.text = element_text(size = 6.5))
pdf(\results/Div-seq/Divseq_results/proportion_region_majorct_time.pdf\, height = 2.7, width = 7.5)
print(plt)
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoibnVsbCBkZXZpY2UgXG4gICAgICAgICAgMSBcbiJ9 -->
null device 1
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Label smoothing
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYWxsX3ByZWRfZiA9IGxpc3QuZmlsZXMoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHBhdHRlcm4gPSBcXHByZWRzX1xcKVxuYWxsX3ByZWRfZiA9IGFsbF9wcmVkX2ZbZ3JlcGwoXFxyZmNcXCwgYWxsX3ByZWRfZikgJiBncmVwbChcXENUXFwsIGFsbF9wcmVkX2YpXVxuIyB0ZXN0IHNldCByZXN1bHRzXG50ZXN0X3JlcyA9IGxhcHBseShhbGxfcHJlZF9mW2dyZXBsKFxcYXguY3N2XFwsIGFsbF9wcmVkX2YpXSwgXG4gICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKHRlc3RfcmVzKSA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQoYWxsX3ByZWRfZltncmVwbChcXGF4LmNzdlxcLCBhbGxfcHJlZF9mKV0sIFxcLlxcLCBmaXhlZCA9IFQpLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSB4WzFdKSlcblxuIyBEaXYtc2VxIHJlc3VsdHNcbmRpdl9yZXMgPSBsYXBwbHkoYWxsX3ByZWRfZltncmVwbChcXGRpdlxcLCBhbGxfcHJlZF9mKV0sIFxuICAgICAgICAgICAgICAgICBmdW5jdGlvbih4KSByZWFkLmNzdihwYXN0ZTAoXFxyZXN1bHRzL0Rpdi1zZXEvXFwsIHgpLCBoZWFkZXIgPSBUKSlcbm5hbWVzKGRpdl9yZXMpID0gdW5saXN0KGxhcHBseShzdHJzcGxpdChhbGxfcHJlZF9mW2dyZXBsKFxcZGl2XFwsIGFsbF9wcmVkX2YpXSwgXFwuXFwsIGZpeGVkID0gVCksXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgeFsxXSkpXG5gYGBcbmBgYCJ9 -->
```r
```r
all_pred_f = list.files(\results/Div-seq/\, pattern = \preds_\)
all_pred_f = all_pred_f[grepl(\rfc\, all_pred_f) & grepl(\CT\, all_pred_f)]
# test set results
test_res = lapply(all_pred_f[grepl(\ax.csv\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(test_res) = unlist(lapply(strsplit(all_pred_f[grepl(\ax.csv\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
# Div-seq results
div_res = lapply(all_pred_f[grepl(\div\, all_pred_f)],
function(x) read.csv(paste0(\results/Div-seq/\, x), header = T))
names(div_res) = unlist(lapply(strsplit(all_pred_f[grepl(\div\, all_pred_f)], \.\, fixed = T),
function(x) x[1]))
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Save new metadata
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuZGl2X2RhdCA9IEZpbmRDbHVzdGVycyhkaXZfZGF0LCByZXNvbHV0aW9uID0gYygyLCA1LCAxMCwgMjAsIDI1LCA1MCksIGFsZ29yaXRobSA9IDIsIHZlcmJvc2UgPSBGKVxuRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxcUk5BX3Nubl9yZXMuMlxcLCBsYWJlbCA9IFQpK05vTGVnZW5kKClcbkRpbVBsb3QoZGl2X2RhdCwgcmVkdWN0aW9uID0gXFx1bWFwX2hhcm1vbnlcXCwgZ3JvdXAuYnkgPSBcXFJOQV9zbm5fcmVzLjVcXCwgbGFiZWwgPSBUKStOb0xlZ2VuZCgpXG5EaW1QbG90KGRpdl9kYXQsIHJlZHVjdGlvbiA9IFxcdW1hcF9oYXJtb255XFwsIGdyb3VwLmJ5ID0gXFxSTkFfc25uX3Jlcy4xMFxcLCBsYWJlbCA9IFQpK05vTGVnZW5kKClcbkRpbVBsb3QoZGl2X2RhdCwgcmVkdWN0aW9uID0gXFx1bWFwX2hhcm1vbnlcXCwgZ3JvdXAuYnkgPSBcXFJOQV9zbm5fcmVzLjIwXFwsIGxhYmVsID0gVCkrTm9MZWdlbmQoKVxuRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxcUk5BX3Nubl9yZXMuMjVcXCwgbGFiZWwgPSBUKStOb0xlZ2VuZCgpXG5EaW1QbG90KGRpdl9kYXQsIHJlZHVjdGlvbiA9IFxcdW1hcF9oYXJtb255XFwsIGdyb3VwLmJ5ID0gXFxSTkFfc25uX3Jlcy41MFxcLCBsYWJlbCA9IFQpK05vTGVnZW5kKClcblxudGFiX292ZXJjbCA9IHRhYmxlKGRpdl9kYXQkUk5BX3Nubl9yZXMuMjAsIGRpdl9kYXQkcHJlZF9jdGFsbClcbnRhYl9vdmVyY2wgPSB0YWJfb3ZlcmNsL3Jvd1N1bXModGFiX292ZXJjbClcbmFwcGx5KHRhYl9vdmVyY2wsIDEsIGZ1bmN0aW9uKHgpIGNvbG5hbWVzKHRhYl9vdmVyY2wpW3doaWNoLm1heCh4KV0pXG5zb3J0KHRhYmxlKGFwcGx5KHRhYl9vdmVyY2wsIDEsIGZ1bmN0aW9uKHgpIGNvbG5hbWVzKHRhYl9vdmVyY2wpW3doaWNoLm1heCh4KV0pKSlcbnNvcnQoYXBwbHkodGFiX292ZXJjbCwgMSwgbWF4KSlcbnN1bShhcHBseSh0YWJfb3ZlcmNsLCAxLCBtYXgpPC41KS9ucm93KHRhYl9vdmVyY2wpXG5cbmRmX292ZXJjbCA9IGRhdGEuZnJhbWUodGFiX292ZXJjbClcbm5vbWFqZGYgPSB1bmlxdWUoZGZfb3ZlcmNsW2RmX292ZXJjbCRGcmVxPC4zMyxjKDEsMSldKVxuY29sbmFtZXMobm9tYWpkZikgPSBjb2xuYW1lcyhkZl9vdmVyY2wpWzE6Ml1cbm1hamRmID0gZGZfb3ZlcmNsW2RmX292ZXJjbCRGcmVxPj0uMzMsMToyXVxubWF0Y2hfZGYgPSByYmluZChtYWpkZiwgbm9tYWpkZlshKG5vbWFqZGYkVmFyMSAlaW4lIG1hamRmJFZhcjEpLF0pXG5cbmRpdl9kYXQkcHJlZF9zbW9vdGhfMjAgPSBtYXRjaF9kZiRWYXIyW21hdGNoKGRpdl9kYXQkUk5BX3Nubl9yZXMuMjAsIG1hdGNoX2RmJFZhcjEpXVxuRGltUGxvdChkaXZfZGF0LCByZWR1Y3Rpb24gPSBcXHVtYXBfaGFybW9ueVxcLCBncm91cC5ieSA9IFxccHJlZF9zbW9vdGhfMjBcXCwgbGFiZWwgPSBUKStOb0xlZ2VuZCgpXG5EaW1QbG90KGRpdl9kYXQsIHJlZHVjdGlvbiA9IFxcdW1hcF9oYXJtb255XFwsIGdyb3VwLmJ5ID0gXFxwcmVkX2N0YWxsXFwsIGxhYmVsID0gVCkrTm9MZWdlbmQoKVxuYGBgXG5gYGAifQ== -->
```r
```r
div_dat = FindClusters(div_dat, resolution = c(2, 5, 10, 20, 25, 50), algorithm = 2, verbose = F)
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.2\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.5\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.10\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.20\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.25\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \RNA_snn_res.50\, label = T)+NoLegend()
tab_overcl = table(div_dat$RNA_snn_res.20, div_dat$pred_ctall)
tab_overcl = tab_overcl/rowSums(tab_overcl)
apply(tab_overcl, 1, function(x) colnames(tab_overcl)[which.max(x)])
sort(table(apply(tab_overcl, 1, function(x) colnames(tab_overcl)[which.max(x)])))
sort(apply(tab_overcl, 1, max))
sum(apply(tab_overcl, 1, max)<.5)/nrow(tab_overcl)
df_overcl = data.frame(tab_overcl)
nomajdf = unique(df_overcl[df_overcl$Freq<.33,c(1,1)])
colnames(nomajdf) = colnames(df_overcl)[1:2]
majdf = df_overcl[df_overcl$Freq>=.33,1:2]
match_df = rbind(majdf, nomajdf[!(nomajdf$Var1 %in% majdf$Var1),])
div_dat$pred_smooth_20 = match_df$Var2[match(div_dat$RNA_snn_res.20, match_df$Var1)]
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \pred_smooth_20\, label = T)+NoLegend()
DimPlot(div_dat, reduction = \umap_harmony\, group.by = \pred_ctall\, label = T)+NoLegend()
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
## Dot plots
Region changes in major predicted cell types
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxubWV0YV9kaXZfcHJlZCA9IGRpdl9kYXRAbWV0YS5kYXRhWyxjKFxcc2FtcGxlXFwsIFxcaGlnaF9sZXZlbF9hbm5vXFwsIFxcaGlnaGxldmVsX2NsdXN0ZXJzXFwsIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFxccHJlZF9yZWdzXFwsIFxccHJlZF9jdGFsbFxcLCBcXHByZWRfc21vb3RoXzIwXFwpXVxud3JpdGUuY3N2KG1ldGFfZGl2X3ByZWQsIGZpbGUgPSBcXHJlc3VsdHMvRGl2LXNlcS9kaXZzZXFfcHJlZGljdGVkX21ldGFkYXRhLmNzdlxcLCBcbiAgICAgICAgICBjb2wubmFtZXMgPSBULCByb3cubmFtZXMgPSBULCBxdW90ZSA9IEYpXG5gYGBcbmBgYCJ9 -->
```r
```r
meta_div_pred = div_dat@meta.data[,c(\sample\, \high_level_anno\, \highlevel_clusters\,
\pred_regs\, \pred_ctall\, \pred_smooth_20\)]
write.csv(meta_div_pred, file = \results/Div-seq/divseq_predicted_metadata.csv\,
col.names = T, row.names = T, quote = F)
<!-- rnb-source-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Scatterpie of previous plot
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxudGFiX2N0X3JlZyA9IHRhYmxlKG1ldGFfZGl2X3ByZWQkc2FtcGxlLFxuICAgICAgICAgICAgICAgICAgIHBhc3RlMChtZXRhX2Rpdl9wcmVkJHByZWRfc21vb3RoXzIwLCBcIi4uXCIsIG1ldGFfZGl2X3ByZWQkcHJlZF9yZWdzKSlcbnRhYl9jdF9yZWcgPSBkYXRhLmZyYW1lKHRhYl9jdF9yZWcpXG5cbnRhYl9jdF9yZWckdGltZSA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQoYXMuY2hhcmFjdGVyKHRhYl9jdF9yZWckVmFyMSksIFwiX1wiKSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHgpIHBhc3RlMCh4WzFdLCBcIl93cGlcIikpKVxudGFiX2N0X3JlZyRwcmVkX3JlZyA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQoYXMuY2hhcmFjdGVyKHRhYl9jdF9yZWckVmFyMiksIFwiLi5cIiwgZml4ZWQgPSBUKSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZ1bmN0aW9uKHgpIHhbMl0pKVxudGFiX2N0X3JlZyRwcmVkX2N0ID0gdW5saXN0KGxhcHBseShzdHJzcGxpdChhcy5jaGFyYWN0ZXIodGFiX2N0X3JlZyRWYXIyKSwgXCIuLlwiLCBmaXhlZCA9IFQpLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgeFsxXSkpXG5cbnRhYl9jdF9yZWckcHJvcCA9IHRhYl9jdF9yZWckRnJlcVxubGFyZ2VzdF9jdCA9IGMoKVxuZm9yKHIgaW4gdW5pcXVlKHRhYl9jdF9yZWckdGltZSkpeyAjIGdldCB0aGUgeCBtb3N0IGZyZXF1ZW50IGNlbGwgdHlwZXMvdHBcbiAgY2MgPSB0YWJfY3RfcmVnJHRpbWU9PXJcbiAgdGFiX2N0X3JlZyRwcm9wW2NjXSA9IHRhYl9jdF9yZWckcHJvcFtjY10vc3VtKHRhYl9jdF9yZWckcHJvcFtjY10pXG4gIGxhcmdlc3RfY3QgPSB1bmlxdWUoYyhsYXJnZXN0X2N0LCBcbiAgICAgICAgICAgICAgICAgICAgICAgIHRhYl9jdF9yZWckcHJlZF9jdFtjY11bb3JkZXIodGFiX2N0X3JlZyRwcm9wW2NjXSwgZGVjcmVhc2luZyA9IFQpXVsxOjhdKSlcbiAgbGFyZ2VzdF9jdCA9IGxhcmdlc3RfY3RbZ3JlcGwoXCJfXCIsIGxhcmdlc3RfY3QpXVxufVxuXG50YWJfY3RfcmVnID0gdGFiX2N0X3JlZ1t0YWJfY3RfcmVnJHByZWRfY3QgJWluJSBsYXJnZXN0X2N0LF1cblxudGFiX3RpbWVjdF9yZWcgPSB0KHJlc2hhcGUyOjpkY2FzdChwcmVkX3JlZ350aW1lK3ByZWRfY3QsIHZhbHVlLnZhciA9IFwiRnJlcVwiLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IHRhYl9jdF9yZWcsIGZpbGwgPSAwKSlcbmNvbG5hbWVzKHRhYl90aW1lY3RfcmVnKSA9IHRhYl90aW1lY3RfcmVnWzEsXVxudGFiX3RpbWVjdF9yZWcgPSBkYXRhLmZyYW1lKHRhYl90aW1lY3RfcmVnWy0xLF0pXG50YWJfdGltZWN0X3JlZyR0aW1lID0gdW5saXN0KGxhcHBseShzdHJzcGxpdChyb3duYW1lcyh0YWJfdGltZWN0X3JlZyksIFwiX3dwaV9cIiksIFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZnVuY3Rpb24oeCkgcGFzdGUwKHhbMV0sIFwiX3dwaVwiKSkpXG50YWJfdGltZWN0X3JlZyRjdCA9IHVubGlzdChsYXBwbHkoc3Ryc3BsaXQocm93bmFtZXModGFiX3RpbWVjdF9yZWcpLCBcIl93cGlfXCIpLCBmdW5jdGlvbih4KSB4WzJdKSlcbnRhYl90aW1lY3RfcmVnJGRvcnNhbCA9IGFzLm51bWVyaWModGFiX3RpbWVjdF9yZWckZG9yc2FsKVxudGFiX3RpbWVjdF9yZWckbGF0ZXJhbCA9IGFzLm51bWVyaWModGFiX3RpbWVjdF9yZWckbGF0ZXJhbClcbnRhYl90aW1lY3RfcmVnJG1lZGlhbCA9IGFzLm51bWVyaWModGFiX3RpbWVjdF9yZWckbWVkaWFsKVxuXG50YWJfdGltZWN0X3JlZyA9IHRhYl90aW1lY3RfcmVnW3RhYl90aW1lY3RfcmVnJGN0IT1cIm1pY3JvZ2xpYV84XCIsXVxuXG5wcm9wc3RpbWUgPSB1bmxpc3QodW5uYW1lKHRhcHBseShhcHBseSh0YWJfdGltZWN0X3JlZ1ssMTozXSwgMSwgc3VtKSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0YWJfdGltZWN0X3JlZyR0aW1lLCBmdW5jdGlvbih4KSB4L3N1bSh4KSkpKVxuXG50YWJfdGltZWN0X3JlZyRwcm9wc3RpbWUgPSBwcm9wc3RpbWVbcm93bmFtZXModGFiX3RpbWVjdF9yZWcpXVxuXG5mb3IociBpbiB1bmlxdWUodGFiX3RpbWVjdF9yZWckY3QpKXtcbiAgY2MgPSB0YWJfdGltZWN0X3JlZyRjdD09clxuICB0YWJfdGltZWN0X3JlZyRwcm9wc3RpbWVbY2NdID0gc2NhbGVzOjpyZXNjYWxlKHRhYl90aW1lY3RfcmVnJHByb3BzdGltZVtjY10sIGMoMCwxKSlcbn1cblxudGFiX3RpbWVjdF9yZWckZ2cgPSAxOm5yb3codGFiX3RpbWVjdF9yZWcpXG5cbnRhYl90aW1lY3RfcmVnJHRpbWUgPSBmYWN0b3IodGFiX3RpbWVjdF9yZWckdGltZSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxldmVscyA9IGMoXCIxX3dwaVwiLFwiMl93cGlcIixcIjRfd3BpXCIsXCI2X3dwaVwiLFwiOF93cGlcIixcIjEyX3dwaVwiKSlcbnRhYl90aW1lY3RfcmVnJGN0ID0gZmFjdG9yKHRhYl90aW1lY3RfcmVnJGN0LCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgIGxldmVscyA9IHJldihjKFwiZXBlbl9jbHVzXzRcIiwgXCJlcGVuX2NsdXNfM1wiLCBcImVwZW5fY2x1c18wXCIsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFwibnBjX1NVQlNFVF83XCIsXCJucGNfU1VCU0VUXzhcIixcIm5wY19TVUJTRVRfM1wiLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBcIm5wY19TVUJTRVRfNFwiLFwibnBjX1NVQlNFVF8xMFwiLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBcImdsdXRfU1VCU0VUXzEwXCIsIFwiZ2x1dF9TVUJTRVRfM1wiLCBcImdsdXRfU1VCU0VUXzFcIixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgXCJnbHV0X1NVQlNFVF80XCIsIFwiZ2x1dF9TVUJTRVRfMFwiLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBcIkdBQkFfU1VCU0VUXzNcIixcIkdBQkFfU1VCU0VUXzRcIixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgXCJHQUJBX1NVQlNFVF8xXCIsXCJHQUJBX1NVQlNFVF8xNFwiKSkpXG5cbnNjcGx0ID0gZ2dwbG90KCkrXG4gIGdlb21fc2NhdHRlcnBpZShhZXMoeD1hcy5udW1lcmljKHRpbWUpLCB5PWFzLm51bWVyaWMoY3QpLCBcbiAgICAgICAgICAgICAgICAgICAgICBncm91cD1nZywgcj1wcm9wc3RpbWUvMiksIGRhdGE9dGFiX3RpbWVjdF9yZWcsXG4gICAgICAgICAgICAgICAgICBjb2xzPWMoXCJtZWRpYWxcIiwgXCJkb3JzYWxcIiwgXCJsYXRlcmFsXCIpKSsgXG4gIHNjYWxlX3hfY29udGludW91cyhsYWJlbHMgPSBsZXZlbHModGFiX3RpbWVjdF9yZWckdGltZSksXG4gICAgICAgICAgICAgICAgICAgICBicmVha3MgPSAxOjYpK1xuICBzY2FsZV95X2NvbnRpbnVvdXMobGFiZWxzID0gbGV2ZWxzKHRhYl90aW1lY3RfcmVnJGN0KSxcbiAgICAgICAgICAgICAgICAgICAgIGJyZWFrcyA9IDE6bGVuZ3RoKGxldmVscyh0YWJfdGltZWN0X3JlZyRjdCkpKStcbiAgc2NhbGVfZmlsbF9tYW51YWwodmFsdWVzID0gcmVnX2NvbHNfc2ltcCkrXG4gIGNvb3JkX2VxdWFsKCkrXG4gIGxhYnMoeCA9IFwidGltZVwiLCB5ID0gXCJjZWxsIHR5cGVcIiwgZmlsbCA9IFwicmVnaW9uXCIpK1xuICB0aGVtZV9idygpK1xuICB0aGVtZShheGlzLnRleHQgPSBlbGVtZW50X3RleHQoY29sb3VyID0gXCJibGFja1wiLCBzaXplID0gNiksXG4gICAgICAgIGF4aXMudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDYuNSksXG4gICAgICAgIGxlZ2VuZC5wb3NpdGlvbiA9IFwibm9uZVwiKVxuXG5wZGYoXCJyZXN1bHRzL0Rpdi1zZXEvRGl2c2VxX3Jlc3VsdHMvcHJvcG9ydGlvbl9yZWdpb25fcHJlZGN0X3RpbWVfcGllLnBkZlwiLCBoZWlnaHQgPSA3LjUsIHdpZHRoID0gMylcbnByaW50KHNjcGx0KVxuZGV2Lm9mZigpXG5gYGAifQ== -->
```r
tab_ct_reg = table(meta_div_pred$sample,
paste0(meta_div_pred$pred_smooth_20, "..", meta_div_pred$pred_regs))
tab_ct_reg = data.frame(tab_ct_reg)
tab_ct_reg$time = unlist(lapply(strsplit(as.character(tab_ct_reg$Var1), "_"),
function(x) paste0(x[1], "_wpi")))
tab_ct_reg$pred_reg = unlist(lapply(strsplit(as.character(tab_ct_reg$Var2), "..", fixed = T),
function(x) x[2]))
tab_ct_reg$pred_ct = unlist(lapply(strsplit(as.character(tab_ct_reg$Var2), "..", fixed = T),
function(x) x[1]))
tab_ct_reg$prop = tab_ct_reg$Freq
largest_ct = c()
for(r in unique(tab_ct_reg$time)){ # get the x most frequent cell types/tp
cc = tab_ct_reg$time==r
tab_ct_reg$prop[cc] = tab_ct_reg$prop[cc]/sum(tab_ct_reg$prop[cc])
largest_ct = unique(c(largest_ct,
tab_ct_reg$pred_ct[cc][order(tab_ct_reg$prop[cc], decreasing = T)][1:8]))
largest_ct = largest_ct[grepl("_", largest_ct)]
}
tab_ct_reg = tab_ct_reg[tab_ct_reg$pred_ct %in% largest_ct,]
tab_timect_reg = t(reshape2::dcast(pred_reg~time+pred_ct, value.var = "Freq",
data = tab_ct_reg, fill = 0))
colnames(tab_timect_reg) = tab_timect_reg[1,]
tab_timect_reg = data.frame(tab_timect_reg[-1,])
tab_timect_reg$time = unlist(lapply(strsplit(rownames(tab_timect_reg), "_wpi_"),
function(x) paste0(x[1], "_wpi")))
tab_timect_reg$ct = unlist(lapply(strsplit(rownames(tab_timect_reg), "_wpi_"), function(x) x[2]))
tab_timect_reg$dorsal = as.numeric(tab_timect_reg$dorsal)
tab_timect_reg$lateral = as.numeric(tab_timect_reg$lateral)
tab_timect_reg$medial = as.numeric(tab_timect_reg$medial)
tab_timect_reg = tab_timect_reg[tab_timect_reg$ct!="microglia_8",]
propstime = unlist(unname(tapply(apply(tab_timect_reg[,1:3], 1, sum),
tab_timect_reg$time, function(x) x/sum(x))))
tab_timect_reg$propstime = propstime[rownames(tab_timect_reg)]
for(r in unique(tab_timect_reg$ct)){
cc = tab_timect_reg$ct==r
tab_timect_reg$propstime[cc] = scales::rescale(tab_timect_reg$propstime[cc], c(0,1))
}
tab_timect_reg$gg = 1:nrow(tab_timect_reg)
tab_timect_reg$time = factor(tab_timect_reg$time,
levels = c("1_wpi","2_wpi","4_wpi","6_wpi","8_wpi","12_wpi"))
tab_timect_reg$ct = factor(tab_timect_reg$ct,
levels = rev(c("epen_clus_4", "epen_clus_3", "epen_clus_0",
"npc_SUBSET_7","npc_SUBSET_8","npc_SUBSET_3",
"npc_SUBSET_4","npc_SUBSET_10",
"glut_SUBSET_10", "glut_SUBSET_3", "glut_SUBSET_1",
"glut_SUBSET_4", "glut_SUBSET_0",
"GABA_SUBSET_3","GABA_SUBSET_4",
"GABA_SUBSET_1","GABA_SUBSET_14")))
scplt = ggplot()+
geom_scatterpie(aes(x=as.numeric(time), y=as.numeric(ct),
group=gg, r=propstime/2), data=tab_timect_reg,
cols=c("medial", "dorsal", "lateral"))+
scale_x_continuous(labels = levels(tab_timect_reg$time),
breaks = 1:6)+
scale_y_continuous(labels = levels(tab_timect_reg$ct),
breaks = 1:length(levels(tab_timect_reg$ct)))+
scale_fill_manual(values = reg_cols_simp)+
coord_equal()+
labs(x = "time", y = "cell type", fill = "region")+
theme_bw()+
theme(axis.text = element_text(colour = "black", size = 6),
axis.title = element_text(size = 6.5),
legend.position = "none")
pdf("results/Div-seq/Divseq_results/proportion_region_predct_time_pie.pdf", height = 7.5, width = 3)
print(scplt)
dev.off()
Scatterpie of major cell types
# count time vs ct..region
major_ct_pred = unlist(lapply(strsplit(meta_div_pred$pred_ctall, "_"), function(x) x[1]))
tab_ct_reg = table(meta_div_pred$sample,
paste0(major_ct_pred, "..", meta_div_pred$pred_regs))
tab_ct_reg = data.frame(tab_ct_reg)
# unravel
tab_ct_reg$time = unlist(lapply(strsplit(as.character(tab_ct_reg$Var1), "_"),
function(x) paste0(x[1], "_wpi")))
tab_ct_reg$pred_reg = unlist(lapply(strsplit(as.character(tab_ct_reg$Var2), "..", fixed = T),
function(x) x[2]))
tab_ct_reg$pred_ct = unlist(lapply(strsplit(as.character(tab_ct_reg$Var2), "..", fixed = T),
function(x) x[1]))
# get proportions
tab_ct_reg$prop = tab_ct_reg$Freq
for(r in unique(tab_ct_reg$time)){
cc = tab_ct_reg$time==r
tab_ct_reg$prop[cc] = tab_ct_reg$prop[cc]/sum(tab_ct_reg$prop[cc])
}
# get counts per region
tab_timect_reg = t(reshape2::dcast(pred_reg~time+pred_ct, value.var = "Freq",
data = tab_ct_reg, fill = 0))
colnames(tab_timect_reg) = tab_timect_reg[1,]
tab_timect_reg = data.frame(tab_timect_reg[-1,])
# unravel
tab_timect_reg$time = unlist(lapply(strsplit(rownames(tab_timect_reg), "_wpi_"),
function(x) paste0(x[1], "_wpi")))
tab_timect_reg$ct = unlist(lapply(strsplit(rownames(tab_timect_reg), "_wpi_"), function(x) x[2]))
tab_timect_reg$dorsal = as.numeric(tab_timect_reg$dorsal)
tab_timect_reg$lateral = as.numeric(tab_timect_reg$lateral)
tab_timect_reg$medial = as.numeric(tab_timect_reg$medial)
# proportions of cell type per time
propstime = unlist(unname(tapply(apply(tab_timect_reg[,1:3], 1, sum),
tab_timect_reg$time, function(x) x/sum(x))))
tab_timect_reg$propstime = propstime[rownames(tab_timect_reg)]
# scale the proportions per cell type
for(r in unique(tab_timect_reg$ct)){
cc = tab_timect_reg$ct==r
tab_timect_reg$propstime[cc] = scales::rescale(tab_timect_reg$propstime[cc], c(0,1))
}
tab_timect_reg$gg = 1:nrow(tab_timect_reg) #required for plot
tab_timect_reg$time = factor(tab_timect_reg$time,
levels = c("1_wpi","2_wpi","4_wpi","6_wpi","8_wpi","12_wpi"))
tab_timect_reg$ct = factor(tab_timect_reg$ct,
levels = rev(c("epen", "npc", "glut", "GABA","oligodendrocyte",
"microglia", "endothelial")))
scplt = ggplot()+
geom_scatterpie(aes(x=as.numeric(time), y=as.numeric(ct),
group=gg, r=propstime/2), data=tab_timect_reg,
cols=c("medial", "dorsal", "lateral"))+
scale_x_continuous(labels = levels(tab_timect_reg$time),
breaks = 1:6)+
scale_y_continuous(labels = levels(tab_timect_reg$ct),
breaks = 1:length(levels(tab_timect_reg$ct)))+
scale_fill_manual(values = reg_cols_simp)+
coord_equal()+
labs(x = "time", y = "cell type", fill = "region")+
theme_bw()+
theme(axis.text = element_text(colour = "black", size = 6),
axis.title = element_text(size = 6.5),
legend.position = "none")
# add global cell type proportions per region
prop_df = data.frame(prop.table(table(major_ct_pred, meta_div_pred$pred_regs)))
prop_df$major_ct_pred = factor(prop_df$major_ct_pred, levels = levels(tab_timect_reg$ct))
col_prop = ggplot()+
geom_col(data = prop_df, mapping = aes(x = Freq, y = major_ct_pred, fill = Var2))+
scale_x_continuous(expand = c(0,0), limits = c(0,.37))+
scale_fill_manual(values = reg_cols_simp)+
labs(x = "proportion")+
theme_void()+
theme(legend.position = "none",
axis.line.y = element_line(),
axis.text.x = element_text(colour = "black", size = 6,
angle = 30, vjust = 1, hjust = 1),
axis.title.x = element_text(size = 6.5),
axis.line.x = element_line(),
axis.ticks.x = element_line())
pdf("results/Div-seq/Divseq_results/proportion_region_predct_time_pie_bar.pdf",
height = 3, width = 3.65)
scplt + col_prop + plot_layout(widths = c(1,0.2)) & theme(plot.margin = unit(c(0,0,0,0), "cm"))
dev.off()
Predicted cell type abundances
tab_ctall = data.frame(table(meta_div_pred$sample,meta_div_pred$pred_ctall))
tab_ctall$majorct = unlist(lapply(strsplit(as.character(tab_ctall$Var2), "_"), function(x) x[1]))
tab_ctall$majorct = factor(tab_ctall$majorct,
levels = c("epen", "npc", "glut", "GABA", "microglia",
"oligodendrocyte", "endothelial"))
ctord = sort(table(meta_div_pred$pred_ctall))
tab_ctall$Var2 = factor(tab_ctall$Var2, levels = names(ctord))
pdf("results/Div-seq/Divseq_results/predCellType_abundance_timepoint.pdf", height = 9.5, width = 7.5)
ggplot(tab_ctall, aes(x = Freq, y = Var2, fill = Var1))+
facet_grid(majorct~Var1, scales = "free", space = "free_y")+
geom_col(colour = "grey50")+
scale_x_continuous(expand = c(0,0))+
scale_fill_manual(values = RColorBrewer::brewer.pal(6, "PuBuGn"))+
theme_bw()+
theme(axis.text.x = element_text(size = 6, colour = "black"),
axis.text.y = element_text(size = 6.5, colour = "black"),
strip.text = element_text(size = 6.5),
legend.position = "none")
dev.off()
pdf("results/Div-seq/Divseq_results/predCellType_abundance_celltype.pdf", height = 9.5, width = 7.5)
ggplot(tab_ctall, aes(x = Freq, y = Var2, fill = Var2))+
facet_grid(majorct~Var1, scales = "free", space = "free_y")+
geom_col()+
scale_x_continuous(expand = c(0,0))+
scale_fill_manual(values = cols_cc)+
theme_bw()+
theme(axis.text.x = element_text(size = 6, colour = "black"),
axis.text.y = element_text(size = 6.5, colour = "black"),
strip.text = element_text(size = 6.5),
legend.position = "none")
dev.off()
Predicted cell type abundances, with predicted region
```r
pdf(\results/Div-seq/Divseq_results/proportion_region_predct_time_pie_bar.pdf\,
height = 3, width = 3.65)
scplt + col_prop + plot_layout(widths = c(1,0.2)) & theme(plot.margin = unit(c(0,0,0,0), \cm\))
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoibnVsbCBkZXZpY2UgXG4gICAgICAgICAgMSBcbiJ9 -->
null device 1
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
Barplots for glut and NPC
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxudGFiX2N0YWxsID0gZGF0YS5mcmFtZSh0YWJsZShtZXRhX2Rpdl9wcmVkJHNhbXBsZSxtZXRhX2Rpdl9wcmVkJHByZWRfY3RhbGwpKVxudGFiX2N0YWxsJG1ham9yY3QgPSB1bmxpc3QobGFwcGx5KHN0cnNwbGl0KGFzLmNoYXJhY3Rlcih0YWJfY3RhbGwkVmFyMiksIFxcX1xcKSwgZnVuY3Rpb24oeCkgeFsxXSkpXG50YWJfY3RhbGwkbWFqb3JjdCA9IGZhY3Rvcih0YWJfY3RhbGwkbWFqb3JjdCwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICBsZXZlbHMgPSBjKFxcZXBlblxcLCBcXG5wY1xcLCBcXGdsdXRcXCwgXFxHQUJBXFwsIFxcbWljcm9nbGlhXFwsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFxcb2xpZ29kZW5kcm9jeXRlXFwsIFxcZW5kb3RoZWxpYWxcXCkpXG5jdG9yZCA9IHNvcnQodGFibGUobWV0YV9kaXZfcHJlZCRwcmVkX2N0YWxsKSlcbnRhYl9jdGFsbCRWYXIyID0gZmFjdG9yKHRhYl9jdGFsbCRWYXIyLCBsZXZlbHMgPSBuYW1lcyhjdG9yZCkpXG5cbnBkZihcXHJlc3VsdHMvRGl2LXNlcS9EaXZzZXFfcmVzdWx0cy9wcmVkQ2VsbFR5cGVfYWJ1bmRhbmNlX3RpbWVwb2ludC5wZGZcXCwgaGVpZ2h0ID0gOS41LCB3aWR0aCA9IDcuNSlcbmdncGxvdCh0YWJfY3RhbGwsIGFlcyh4ID0gRnJlcSwgeSA9IFZhcjIsIGZpbGwgPSBWYXIxKSkrXG4gIGZhY2V0X2dyaWQobWFqb3JjdH5WYXIxLCBzY2FsZXMgPSBcXGZyZWVcXCwgc3BhY2UgPSBcXGZyZWVfeVxcKStcbiAgZ2VvbV9jb2woY29sb3VyID0gXFxncmV5NTBcXCkrXG4gIHNjYWxlX3hfY29udGludW91cyhleHBhbmQgPSBjKDAsMCkpK1xuICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBSQ29sb3JCcmV3ZXI6OmJyZXdlci5wYWwoNiwgXFxQdUJ1R25cXCkpK1xuICB0aGVtZV9idygpK1xuICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChzaXplID0gNiwgY29sb3VyID0gXFxibGFja1xcKSxcbiAgICAgICAgYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDYuNSwgY29sb3VyID0gXFxibGFja1xcKSxcbiAgICAgICAgc3RyaXAudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNi41KSxcbiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gXFxub25lXFwpXG5kZXYub2ZmKClcbmBgYFxuYGBgIn0= -->
```r
```r
tab_ctall = data.frame(table(meta_div_pred$sample,meta_div_pred$pred_ctall))
tab_ctall$majorct = unlist(lapply(strsplit(as.character(tab_ctall$Var2), \_\), function(x) x[1]))
tab_ctall$majorct = factor(tab_ctall$majorct,
levels = c(\epen\, \npc\, \glut\, \GABA\, \microglia\,
\oligodendrocyte\, \endothelial\))
ctord = sort(table(meta_div_pred$pred_ctall))
tab_ctall$Var2 = factor(tab_ctall$Var2, levels = names(ctord))
pdf(\results/Div-seq/Divseq_results/predCellType_abundance_timepoint.pdf\, height = 9.5, width = 7.5)
ggplot(tab_ctall, aes(x = Freq, y = Var2, fill = Var1))+
facet_grid(majorct~Var1, scales = \free\, space = \free_y\)+
geom_col(colour = \grey50\)+
scale_x_continuous(expand = c(0,0))+
scale_fill_manual(values = RColorBrewer::brewer.pal(6, \PuBuGn\))+
theme_bw()+
theme(axis.text.x = element_text(size = 6, colour = \black\),
axis.text.y = element_text(size = 6.5, colour = \black\),
strip.text = element_text(size = 6.5),
legend.position = \none\)
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoibnVsbCBkZXZpY2UgXG4gICAgICAgICAgMSBcbiJ9 -->
null device 1
<!-- rnb-output-end -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucGRmKFxccmVzdWx0cy9EaXYtc2VxL0RpdnNlcV9yZXN1bHRzL3ByZWRDZWxsVHlwZV9hYnVuZGFuY2VfY2VsbHR5cGUucGRmXFwsIGhlaWdodCA9IDkuNSwgd2lkdGggPSA3LjUpXG5gYGBcbmBgYHJcbmdncGxvdCh0YWJfY3RhbGwsIGFlcyh4ID0gRnJlcSwgeSA9IFZhcjIsIGZpbGwgPSBWYXIyKSkrXG4gIGZhY2V0X2dyaWQobWFqb3JjdH5WYXIxLCBzY2FsZXMgPSBcXGZyZWVcXCwgc3BhY2UgPSBcXGZyZWVfeVxcKStcbiAgZ2VvbV9jb2woKStcbiAgc2NhbGVfeF9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwwKSkrXG4gIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbHNfY2MpK1xuICB0aGVtZV9idygpK1xuICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChzaXplID0gNiwgY29sb3VyID0gXFxibGFja1xcKSxcbiAgICAgICAgYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDYuNSwgY29sb3VyID0gXFxibGFja1xcKSxcbiAgICAgICAgc3RyaXAudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNi41KSxcbiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gXFxub25lXFwpXG5kZXYub2ZmKClcbmBgYFxuYGBgIn0= -->
```r
```r
pdf(\results/Div-seq/Divseq_results/predCellType_abundance_celltype.pdf\, height = 9.5, width = 7.5)
ggplot(tab_ctall, aes(x = Freq, y = Var2, fill = Var2))+
facet_grid(majorct~Var1, scales = \free\, space = \free_y\)+
geom_col()+
scale_x_continuous(expand = c(0,0))+
scale_fill_manual(values = cols_cc)+
theme_bw()+
theme(axis.text.x = element_text(size = 6, colour = \black\),
axis.text.y = element_text(size = 6.5, colour = \black\),
strip.text = element_text(size = 6.5),
legend.position = \none\)
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoibnVsbCBkZXZpY2UgXG4gICAgICAgICAgMSBcbiJ9 -->
null device 1
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
## UMAPs
UMAP with predicted cell types
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucGxvdF9kZiA9IG1ldGFfZGl2X3ByZWRbLGMoMSw0LDUpXVxucGxvdF9kZiRzYW1wbGVbcGxvdF9kZiRzYW1wbGU9PVwiMV93cGlfcG9zXCJdID0gXCIxX3dwaVwiXG5wbG90X2RmJHNhbXBsZSA9IGZhY3RvcihwbG90X2RmJHNhbXBsZSwgbGV2ZWxzID0gcmV2KHBhc3RlMChjKDEsMiw0LDYsOCwxMiksIFwiX3dwaVwiKSkpXG5wbG90X2RmJHByZWRfcmVncyA9IGZhY3RvcihwbG90X2RmJHByZWRfcmVncywgbGV2ZWxzID0gYyhcIm1lZGlhbFwiLCBcImRvcnNhbFwiLCBcImxhdGVyYWxcIikpXG5mb3IoaSBpbiBjKFwiZ2x1dF9cIiwgXCJucGNfXCIpKXtcbiAgc3ViX3Bsb3RfZGYgPSBwbG90X2RmW2dyZXBsKGksIHBsb3RfZGYkcHJlZF9jdGFsbCksXVxuICBzdWJfcGxvdF9kZiRwcmVkX2N0YWxsID0gZ3N1YihwYXN0ZTAoaSwgXCJTVUJTRVRfXCIpLCBcIlwiLCBzdWJfcGxvdF9kZiRwcmVkX2N0YWxsKVxuICBsdmxzID0gaWYoaT09XCJnbHV0X1wiKXtcbiAgICBjKDksMTAsMjIsMjcsMTksMTQsMTcsMjYsOCw2LDIwLDEsMTEsMTMsMjUsMiwyMSwxNiwzLDAsMTIsNyw1LDQsMjQsMTUsMTgpXG4gIH0gZWxzZSB7YygyLDEzLDAsMTQsMSwzLDksNCwxMSw3KX1cbiAgc3ViX3Bsb3RfZGYgPSBzdWJfcGxvdF9kZltzdWJfcGxvdF9kZiRwcmVkX2N0YWxsICVpbiUgYXMuY2hhcmFjdGVyKGx2bHMpLF1cbiAgc3ViX3Bsb3RfZGYkcHJlZF9jdGFsbCA9IGZhY3RvcihzdWJfcGxvdF9kZiRwcmVkX2N0YWxsLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBsZXZlbHMgPSBhcy5jaGFyYWN0ZXIobHZscykpXG4gIGxpbW1heCA9IG1heCh0YWJsZShzdWJfcGxvdF9kZiRwcmVkX2N0YWxsKSlcbiAgXG4gIGZvcihnZyBpbiBjKFwic2FtcGxlXCIsIFwicHJlZF9yZWdzXCIsIFwiZ3JleTI1XCIpKXtcbiAgICBjb2xzX3VzZSA9IGlmKGdnPT1cInByZWRfcmVnc1wiKXtcbiAgICAgIHJlZ19jb2xzX3NpbXBcbiAgICB9IGVsc2UgaWYoZ2c9PVwic2FtcGxlXCIpe1xuICAgICAgcmV2KGNvbG9yUmFtcFBhbGV0dGUoYyhcImdyZXk4NVwiLCBcIiNkZjE2ZGZcIikpKDYpKVxuICAgIH1cbiAgICBwbHQgPSBnZ3Bsb3Qoc3ViX3Bsb3RfZGYpK1xuICAgICAgc2NhbGVfeV9jb250aW51b3VzKGV4cGFuZCA9IGMoMCwwKSwgbGltaXRzID0gYygwLCBsaW1tYXgrMC4wNSpsaW1tYXgpKStcbiAgICAgIHRoZW1lX2NsYXNzaWMoKStcbiAgICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGNvbG91ciA9IFwiYmxhY2tcIiwgc2l6ZSA9IDcpLFxuICAgICAgICAgICAgYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoY29sb3VyID0gXCJibGFja1wiLCBzaXplID0gNiksXG4gICAgICAgICAgICBheGlzLnRpdGxlLnggPSBlbGVtZW50X2JsYW5rKCksXG4gICAgICAgICAgICBheGlzLnRpdGxlLnkgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDcpLFxuICAgICAgICAgICAgbGVnZW5kLmtleS5zaXplID0gdW5pdCgwLjQsIFwiY21cIiksXG4gICAgICAgICAgICBsZWdlbmQudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNiksXG4gICAgICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDcpKVxuICAgIGlmKGdnIT1cImdyZXkyNVwiKXtcbiAgICAgIHBsdCA9IHBsdCtnZW9tX2JhcihtYXBwaW5nID0gYWVzX3N0cmluZyh4ID0gXCJwcmVkX2N0YWxsXCIsIGZpbGwgPSBnZykpK1xuICAgICAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2xzX3VzZSlcbiAgICB9IGVsc2V7XG4gICAgICBwbHQgPSBwbHQrZ2VvbV9iYXIobWFwcGluZyA9IGFlc19zdHJpbmcoeCA9IFwicHJlZF9jdGFsbFwiKSwgZmlsbCA9IGdnKVxuICAgIH1cbiAgICBwZGYocGFzdGUwKFwicmVzdWx0cy9EaXYtc2VxL2JhcnBsb3RzX2dsdXRucGMvYmFyX1wiLCBpLCBnZywgXCIucGRmXCIpLCB3aWR0aCA9IDQuOCwgaGVpZ2h0ID0gMS4yNSlcbiAgICBwcmludChwbHQpXG4gICAgZGV2Lm9mZigpXG4gIH1cbn1cbmBgYCJ9 -->
```r
plot_df = meta_div_pred[,c(1,4,5)]
plot_df$sample[plot_df$sample=="1_wpi_pos"] = "1_wpi"
plot_df$sample = factor(plot_df$sample, levels = rev(paste0(c(1,2,4,6,8,12), "_wpi")))
plot_df$pred_regs = factor(plot_df$pred_regs, levels = c("medial", "dorsal", "lateral"))
for(i in c("glut_", "npc_")){
sub_plot_df = plot_df[grepl(i, plot_df$pred_ctall),]
sub_plot_df$pred_ctall = gsub(paste0(i, "SUBSET_"), "", sub_plot_df$pred_ctall)
lvls = if(i=="glut_"){
c(9,10,22,27,19,14,17,26,8,6,20,1,11,13,25,2,21,16,3,0,12,7,5,4,24,15,18)
} else {c(2,13,0,14,1,3,9,4,11,7)}
sub_plot_df = sub_plot_df[sub_plot_df$pred_ctall %in% as.character(lvls),]
sub_plot_df$pred_ctall = factor(sub_plot_df$pred_ctall,
levels = as.character(lvls))
limmax = max(table(sub_plot_df$pred_ctall))
for(gg in c("sample", "pred_regs", "grey25")){
cols_use = if(gg=="pred_regs"){
reg_cols_simp
} else if(gg=="sample"){
rev(colorRampPalette(c("grey85", "#df16df"))(6))
}
plt = ggplot(sub_plot_df)+
scale_y_continuous(expand = c(0,0), limits = c(0, limmax+0.05*limmax))+
theme_classic()+
theme(axis.text.x = element_text(colour = "black", size = 7),
axis.text.y = element_text(colour = "black", size = 6),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 7),
legend.key.size = unit(0.4, "cm"),
legend.text = element_text(size = 6),
legend.title = element_text(size = 7))
if(gg!="grey25"){
plt = plt+geom_bar(mapping = aes_string(x = "pred_ctall", fill = gg))+
scale_fill_manual(values = cols_use)
} else{
plt = plt+geom_bar(mapping = aes_string(x = "pred_ctall"), fill = gg)
}
pdf(paste0("results/Div-seq/barplots_glutnpc/bar_", i, gg, ".pdf"), width = 4.8, height = 1.25)
print(plt)
dev.off()
}
}
Line plot for major cell types
```r
tab_ctall = data.frame(table(meta_div_pred$sample,
paste0(meta_div_pred$pred_ctall, \..\, meta_div_pred$pred_regs)))
tab_ctall$ct = unlist(lapply(strsplit(as.character(tab_ctall$Var2), \..\, fixed = T),
function(x) x[1]))
tab_ctall$reg = unlist(lapply(strsplit(as.character(tab_ctall$Var2), \..\, fixed = T),
function(x) x[2]))
tab_ctall$majorct = unlist(lapply(strsplit(as.character(tab_ctall$Var2), \_\), function(x) x[1]))
tab_ctall$majorct = factor(tab_ctall$majorct,
levels = c(\epen\, \npc\, \glut\, \GABA\, \microglia\,
\oligodendrocyte\, \endothelial\))
ctord = sort(table(meta_div_pred$pred_ctall))
tab_ctall$ct = factor(tab_ctall$ct, levels = names(ctord))
pdf(\results/Div-seq/Divseq_results/predCellType_abundance_predRegion.pdf\, height = 9.5, width = 7.5)
ggplot(tab_ctall, aes(x = Freq, y = ct, fill = reg))+
facet_grid(majorct~Var1, scales = \free\, space = \free_y\)+
geom_col()+
scale_x_continuous(expand = c(0,0))+
scale_fill_manual(values = reg_cols_simp)+
theme_bw()+
theme(axis.text.x = element_text(size = 6, colour = \black\),
axis.text.y = element_text(size = 6.5, colour = \black\),
strip.text = element_text(size = 6.5),
legend.position = \none\)
dev.off()
<!-- rnb-source-end -->
<!-- rnb-output-begin eyJkYXRhIjoibnVsbCBkZXZpY2UgXG4gICAgICAgICAgMSBcbiJ9 -->
null device 1
<!-- rnb-output-end -->
<!-- rnb-chunk-end -->
<!-- rnb-text-begin -->
# Integrating SS and 1 wpi neg
Load data
<!-- rnb-text-end -->
<!-- rnb-chunk-begin -->
<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxubWFqb3JjdCA9IHNhcHBseShzdHJzcGxpdChtZXRhX2Rpdl9wcmVkJHByZWRfY3RhbGwsIFwiX1wiKSwgZnVuY3Rpb24oeCkgeFsxXSlcbm1ham9yY3RbbWV0YV9kaXZfcHJlZCRwcmVkX2N0YWxsICVpbiUgYyhcImVwZW5fY2x1c180XCIsIFwiZXBlbl9jbHVzXzNcIildID0gXCJlcGVuX2FcIlxubWFqb3JjdFttZXRhX2Rpdl9wcmVkJHByZWRfY3RhbGwgJWluJSBjKFwiZXBlbl9jbHVzXzFcIiwgXCJlcGVuX2NsdXNfMTNcIiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgXCJlcGVuX2NsdXNfN1wiLFwiZXBlbl9jbHVzXzE0XCIpXSA9IFwiZXBlbl9uXCJcbnBsb3RfZGYgPSByZXNoYXBlMjo6bWVsdChwcm9wLnRhYmxlKHRhYmxlKG1ham9yY3QsIG1ldGFfZGl2X3ByZWQkc2FtcGxlKSwgbWFyZ2luID0gMikpXG5wbG90X2RmJFZhcjIgPSBhcy5jaGFyYWN0ZXIocGxvdF9kZiRWYXIyKVxucGxvdF9kZiRWYXIyW3Bsb3RfZGYkVmFyMj09XCIxX3dwaV9wb3NcIl0gPSBcIjFfd3BpXCJcbnBsb3RfZGYkVmFyMiA9IGZhY3RvcihwbG90X2RmJFZhcjIsIGxldmVscyA9IGMoXCIxX3dwaVwiLCBcIjJfd3BpXCIsIFwiNF93cGlcIiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIFwiNl93cGlcIiwgXCI4X3dwaVwiLCBcIjEyX3dwaVwiKSlcblxuY29sc191c2UgPSBjKHVubmFtZShjb2xzX2NjW2dyZXBsKFwiZXBlblwiLCBuYW1lcyhjb2xzX2NjKSldW2MoMSwxNSw4KV0pLCBcIiM2MDkxYmFcIiwgXCIjZDA1OTU5XCIsXG4gICAgICAgICAgICAgXCIjZmViNzJhXCIsIFwiI0U2NTMwRFwiLCBcIiNFNDNEODhcIiwgXCIjNzEyMTY2XCIpXG5uYW1lcyhjb2xzX3VzZSkgPSBjKFwiZXBlbl9hXCIsIFwiZXBlbl9uXCIsIFwiZXBlbl9xXCIsIFwiZ2x1dFwiLCBcIkdBQkFcIiwgXCJucGNcIiwgXG4gICAgICAgICAgICAgICAgICAgIFwibWljcm9nbGlhXCIsIFwib2xpZ29kZW5kcm9jeXRlXCIsIFwiZW5kb3RoZWxpYWxcIilcbmxuX3BsdCA9IGdncGxvdChwbG90X2RmLCBhZXMoeCA9IFZhcjIsIHkgPSB2YWx1ZSwgXG4gICAgICAgICAgICAgICAgICAgIGdyb3VwID0gbWFqb3JjdCwgY29sb3VyID0gbWFqb3JjdCkpK1xuICBnZW9tX3BvaW50KHNpemUgPSAxKStcbiAgZ2VvbV9saW5lKCkrXG4gIHNjYWxlX2NvbG91cl9tYW51YWwodmFsdWVzID0gY29sc191c2UpK1xuICBsYWJzKHggPSBcInRpbWUgKHdwaSlcIiwgeSA9IFwicHJvcG9ydGlvbiBwZXIgdGltZXBvaW50XCIsXG4gICAgICAgY29sb3VyID0gXCJtYWpvciBjZWxsIHR5cGVcXG4ocHJlZGljdGVkKVwiKStcbiAgdGhlbWVfY2xhc3NpYygpK1xuICB0aGVtZShhc3BlY3QucmF0aW8gPSAxLzIsXG4gICAgICAgIGF4aXMudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplID0gNywgY29sb3VyID0gXCJibGFja1wiKSxcbiAgICAgICAgYXhpcy50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplID0gNy41KSxcbiAgICAgICAgbGVnZW5kLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDcpLFxuICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDcuNSksXG4gICAgICAgIGxlZ2VuZC5rZXkuc2l6ZSA9IHVuaXQoMC41NSwgXCJjbVwiKSlcblxucGRmKFwicmVzdWx0cy9EaXYtc2VxL0RpdnNlcV9yZXN1bHRzL2RpdnNlcV9tYWpvcmN0X2xpbmVzLnBkZlwiLCBoZWlnaHQgPSAzLCB3aWR0aCA9IDUpXG5wcmludChsbl9wbHQpXG5kZXYub2ZmKClcbmBgYCJ9 -->
```r
majorct = sapply(strsplit(meta_div_pred$pred_ctall, "_"), function(x) x[1])
majorct[meta_div_pred$pred_ctall %in% c("epen_clus_4", "epen_clus_3")] = "epen_a"
majorct[meta_div_pred$pred_ctall %in% c("epen_clus_1", "epen_clus_13",
"epen_clus_7","epen_clus_14")] = "epen_n"
plot_df = reshape2::melt(prop.table(table(majorct, meta_div_pred$sample), margin = 2))
plot_df$Var2 = as.character(plot_df$Var2)
plot_df$Var2[plot_df$Var2=="1_wpi_pos"] = "1_wpi"
plot_df$Var2 = factor(plot_df$Var2, levels = c("1_wpi", "2_wpi", "4_wpi",
"6_wpi", "8_wpi", "12_wpi"))
cols_use = c(unname(cols_cc[grepl("epen", names(cols_cc))][c(1,15,8)]), "#6091ba", "#d05959",
"#feb72a", "#E6530D", "#E43D88", "#712166")
names(cols_use) = c("epen_a", "epen_n", "epen_q", "glut", "GABA", "npc",
"microglia", "oligodendrocyte", "endothelial")
ln_plt = ggplot(plot_df, aes(x = Var2, y = value,
group = majorct, colour = majorct))+
geom_point(size = 1)+
geom_line()+
scale_colour_manual(values = cols_use)+
labs(x = "time (wpi)", y = "proportion per timepoint",
colour = "major cell type\n(predicted)")+
theme_classic()+
theme(aspect.ratio = 1/2,
axis.text = element_text(size = 7, colour = "black"),
axis.title = element_text(size = 7.5),
legend.text = element_text(size = 7),
legend.title = element_text(size = 7.5),
legend.key.size = unit(0.55, "cm"))
pdf("results/Div-seq/Divseq_results/divseq_majorct_lines.pdf", height = 3, width = 5)
print(ln_plt)
dev.off()
null device
1
Merge datasets
ax_srat = readRDS("data/expression/axolotl_reclust/all_nuclei_clustered_highlevel_anno.RDS")
meta = read.csv("data/annotations/axolotl_all_umeta.csv",
header = T, row.names = 1)
ax_srat = AddMetaData(ax_srat, metadata = meta)
meta_regs = read.csv("data/processed/multiome/WP_region_predictions.csv",
header = T, row.names = 1)
newcellnames = rownames(meta_regs)
newcellnames = gsub("-a1-1", "-1_1", newcellnames)
newcellnames = gsub("-a1-2", "-1_2", newcellnames)
newcellnames = gsub("-a3-1", "-1_3", newcellnames)
newcellnames = gsub("-a3-2", "-1_4", newcellnames)
rownames(meta_regs) = newcellnames
ax_srat = AddMetaData(ax_srat, metadata = meta_regs)
ax_srat$regions_all = ax_srat$pred_regions_top
ax_srat$regions_all[is.na(ax_srat$regions_all)] = ax_srat$regions[is.na(ax_srat$regions_all)]
neg_srat = readRDS("data/expression/axolotl_reclust/divseq_1wpi_neg.RDS")
reg_pred = read.csv("results/Div-seq/preds_lr_regions_neg_all.csv", header = T, row.names = 1)
ct_pred = read.csv("results/Div-seq/preds_rfc_CT_neg_all.csv", header = T, row.names = 1)
neg_srat$pred_regs = reg_pred[colnames(neg_srat),"pred_lr"]
neg_srat$pred_ctall = ct_pred[colnames(neg_srat),"pred_rfc"]
Joint processing
# filter metadata
ax_srat@meta.data = ax_srat@meta.data[,c("nCount_RNA", "nFeature_RNA", "percent.mt", "chem", "high_level_anno", "cellclusters", "regions_all")]
colnames(ax_srat@meta.data)[4] = "sample"
neg_srat@meta.data = neg_srat@meta.data[,c("nCount_RNA", "nFeature_RNA", "percent.mt", "sample", "seurat_clusters", "pred_ctall", "pred_regs"), drop = F]
ax_neg_srat = merge(ax_srat, neg_srat)
Run UMAP
ax_neg_srat$sample_dat = ax_neg_srat$sample
ax_neg_srat$sample_dat[grepl("wpi", ax_neg_srat$sample_dat)] = "div"
ax_neg_srat = NormalizeData(ax_neg_srat)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
ax_neg_srat = FindVariableFeatures(ax_neg_srat, nfeatures = 10000)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
ax_neg_srat = ScaleData(ax_neg_srat, vars.to.regress = c("nCount_RNA", "percent.mt", "sample_dat"))
Regressing out nCount_RNA, percent.mt, sample_dat
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 3%
|
|=== | 4%
|
|==== | 4%
|
|==== | 5%
|
|===== | 5%
|
|===== | 6%
|
|====== | 6%
|
|====== | 7%
|
|======= | 7%
|
|======= | 8%
|
|======= | 9%
|
|======== | 9%
|
|======== | 10%
|
|========= | 10%
|
|========= | 11%
|
|========== | 11%
|
|========== | 12%
|
|=========== | 12%
|
|=========== | 13%
|
|============ | 13%
|
|============ | 14%
|
|============= | 14%
|
|============= | 15%
|
|============= | 16%
|
|============== | 16%
|
|============== | 17%
|
|=============== | 17%
|
|=============== | 18%
|
|================ | 18%
|
|================ | 19%
|
|================= | 19%
|
|================= | 20%
|
|================== | 20%
|
|================== | 21%
|
|=================== | 21%
|
|=================== | 22%
|
|==================== | 22%
|
|==================== | 23%
|
|==================== | 24%
|
|===================== | 24%
|
|===================== | 25%
|
|====================== | 25%
|
|====================== | 26%
|
|======================= | 26%
|
|======================= | 27%
|
|======================== | 27%
|
|======================== | 28%
|
|========================= | 28%
|
|========================= | 29%
|
|========================== | 29%
|
|========================== | 30%
|
|=========================== | 30%
|
|=========================== | 31%
|
|=========================== | 32%
|
|============================ | 32%
|
|============================ | 33%
|
|============================= | 33%
|
|============================= | 34%
|
|============================== | 34%
|
|============================== | 35%
|
|=============================== | 35%
|
|=============================== | 36%
|
|================================ | 36%
|
|================================ | 37%
|
|================================= | 37%
|
|================================= | 38%
|
|================================== | 39%
|
|================================== | 40%
|
|=================================== | 40%
|
|=================================== | 41%
|
|==================================== | 41%
|
|==================================== | 42%
|
|===================================== | 42%
|
|===================================== | 43%
|
|====================================== | 43%
|
|====================================== | 44%
|
|======================================= | 44%
|
|======================================= | 45%
|
|======================================== | 45%
|
|======================================== | 46%
|
|======================================== | 47%
|
|========================================= | 47%
|
|========================================= | 48%
|
|========================================== | 48%
|
|========================================== | 49%
|
|=========================================== | 49%
|
|=========================================== | 50%
|
|============================================ | 50%
|
|============================================ | 51%
|
|============================================= | 51%
|
|============================================= | 52%
|
|============================================== | 52%
|
|============================================== | 53%
|
|=============================================== | 53%
|
|=============================================== | 54%
|
|=============================================== | 55%
|
|================================================ | 55%
|
|================================================ | 56%
|
|================================================= | 56%
|
|================================================= | 57%
|
|================================================== | 57%
|
|================================================== | 58%
|
|=================================================== | 58%
|
|=================================================== | 59%
|
|==================================================== | 59%
|
|==================================================== | 60%
|
|===================================================== | 60%
|
|===================================================== | 61%
|
|====================================================== | 62%
|
|====================================================== | 63%
|
|======================================================= | 63%
|
|======================================================= | 64%
|
|======================================================== | 64%
|
|======================================================== | 65%
|
|========================================================= | 65%
|
|========================================================= | 66%
|
|========================================================== | 66%
|
|========================================================== | 67%
|
|=========================================================== | 67%
|
|=========================================================== | 68%
|
|============================================================ | 68%
|
|============================================================ | 69%
|
|============================================================ | 70%
|
|============================================================= | 70%
|
|============================================================= | 71%
|
|============================================================== | 71%
|
|============================================================== | 72%
|
|=============================================================== | 72%
|
|=============================================================== | 73%
|
|================================================================ | 73%
|
|================================================================ | 74%
|
|================================================================= | 74%
|
|================================================================= | 75%
|
|================================================================== | 75%
|
|================================================================== | 76%
|
|=================================================================== | 76%
|
|=================================================================== | 77%
|
|=================================================================== | 78%
|
|==================================================================== | 78%
|
|==================================================================== | 79%
|
|===================================================================== | 79%
|
|===================================================================== | 80%
|
|====================================================================== | 80%
|
|====================================================================== | 81%
|
|======================================================================= | 81%
|
|======================================================================= | 82%
|
|======================================================================== | 82%
|
|======================================================================== | 83%
|
|========================================================================= | 83%
|
|========================================================================= | 84%
|
|========================================================================== | 84%
|
|========================================================================== | 85%
|
|========================================================================== | 86%
|
|=========================================================================== | 86%
|
|=========================================================================== | 87%
|
|============================================================================ | 87%
|
|============================================================================ | 88%
|
|============================================================================= | 88%
|
|============================================================================= | 89%
|
|============================================================================== | 89%
|
|============================================================================== | 90%
|
|=============================================================================== | 90%
|
|=============================================================================== | 91%
|
|================================================================================ | 91%
|
|================================================================================ | 92%
|
|================================================================================ | 93%
|
|================================================================================= | 93%
|
|================================================================================= | 94%
|
|================================================================================== | 94%
|
|================================================================================== | 95%
|
|=================================================================================== | 95%
|
|=================================================================================== | 96%
|
|==================================================================================== | 96%
|
|==================================================================================== | 97%
|
|===================================================================================== | 97%
|
|===================================================================================== | 98%
|
|====================================================================================== | 98%
|
|====================================================================================== | 99%
|
|=======================================================================================| 99%
|
|=======================================================================================| 100%
Centering and scaling data matrix
|
| | 0%
|
|========= | 10%
|
|================= | 20%
|
|========================== | 30%
|
|=================================== | 40%
|
|============================================ | 50%
|
|==================================================== | 60%
|
|============================================================= | 70%
|
|====================================================================== | 80%
|
|============================================================================== | 90%
|
|=======================================================================================| 100%
ax_neg_srat = RunPCA(ax_neg_srat, verbose = F)
ElbowPlot(ax_neg_srat, ndims = 50)
Integrate with Harmony based on multiome/v3.1/div-seq
ax_neg_srat = RunUMAP(ax_neg_srat, dims = 1:20, verbose = F)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
UMAPPlot(ax_neg_srat, group.by = "sample_dat")
UMAPPlot(ax_neg_srat, group.by = "cellclusters")
Check integration
ax_neg_srat$sample_simp = ax_neg_srat$sample
ax_neg_srat$sample_simp[ax_neg_srat$sample_simp!="1_wpi_neg"] = "SS"
ax_neg_srat$ctlabels = ax_neg_srat$cellclusters
ax_neg_srat$ctlabels[is.na(ax_neg_srat$ctlabels)] = ax_neg_srat$pred_ctall[is.na(ax_neg_srat$ctlabels)]
ax_neg_srat$reglabels = ax_neg_srat$regions_all
ax_neg_srat$reglabels[is.na(ax_neg_srat$reglabels)] = ax_neg_srat$pred_regs[is.na(ax_neg_srat$reglabels)]
DimPlot(ax_neg_srat, group.by = "sample_dat", reduction = "umap_harmony")
DimPlot(ax_neg_srat, group.by = "ctlabels", reduction = "umap_harmony")
plt1 = DimPlot(ax_neg_srat, group.by = "ctlabels", reduction = "umap_harmony",
label = T, split.by = "sample_dat")
plt2 = DimPlot(ax_neg_srat, group.by = "ctlabels", reduction = "umap_harmony",
label = T, split.by = "sample")
plt3 = DimPlot(ax_neg_srat, group.by = "sample_dat", reduction = "umap_harmony", label = T)
plt4 = DimPlot(ax_neg_srat, group.by = "reglabels", reduction = "umap_harmony", label = T)
plt5 = DimPlot(ax_neg_srat, group.by = "ctlabels", reduction = "umap_harmony",
label = T, split.by = "reglabels")
DimPlot(ax_neg_srat, group.by = "ctlabels", reduction = "umap_harmony",
label = T, split.by = "sample_simp")+NoLegend()
DimPlot(ax_neg_srat, group.by = "sample_simp", reduction = "umap_harmony")
Clustering
DE between protocols
Filter the DE genes to those significant and meaningful
ct_mk_l_filt = list()
for(n in names(ct_mk_l)){
ct_mk_l_filt[[n]] = ct_mk_l[[n]][ct_mk_l[[n]]$padj<=0.05 & ct_mk_l[[n]]$logFC>=0.3,]
ct_mk_l_filt[[n]] = ct_mk_l_filt[[n]][!startsWith(ct_mk_l_filt[[n]]$feature, "AMEX") &
!startsWith(ct_mk_l_filt[[n]]$feature, "LOC") &
!grepl("..", ct_mk_l_filt[[n]]$feature, fixed = T),]
}
cl_mk_l_filt = list()
for(n in names(cl_mk_l)){
cl_mk_l_filt[[n]] = cl_mk_l[[n]][cl_mk_l[[n]]$padj<=0.05 & cl_mk_l[[n]]$logFC>=0.3,]
cl_mk_l_filt[[n]] = cl_mk_l_filt[[n]][!startsWith(cl_mk_l_filt[[n]]$feature, "AMEX") &
!startsWith(cl_mk_l_filt[[n]]$feature, "LOC") &
!grepl("..", cl_mk_l_filt[[n]]$feature, fixed = T),]
}
Count DE genes per cell type and condition
ct_mk_l_filt = list()
for(n in names(ct_mk_l)){
ct_mk_l_filt[[n]] = ct_mk_l[[n]][ct_mk_l[[n]]$padj<=0.05 & ct_mk_l[[n]]$logFC>=0.3,]
ct_mk_l_filt[[n]] = ct_mk_l_filt[[n]][!startsWith(ct_mk_l_filt[[n]]$feature, "AMEX") &
!startsWith(ct_mk_l_filt[[n]]$feature, "LOC") &
!grepl("..", ct_mk_l_filt[[n]]$feature, fixed = T),]
}
cl_mk_l_filt = list()
for(n in names(cl_mk_l)){
cl_mk_l_filt[[n]] = cl_mk_l[[n]][cl_mk_l[[n]]$padj<=0.05 & cl_mk_l[[n]]$logFC>=0.3,]
cl_mk_l_filt[[n]] = cl_mk_l_filt[[n]][!startsWith(cl_mk_l_filt[[n]]$feature, "AMEX") &
!startsWith(cl_mk_l_filt[[n]]$feature, "LOC") &
!grepl("..", cl_mk_l_filt[[n]]$feature, fixed = T),]
}
Get GO Terms, without repeated genes
n_de = matrix(0, ncol = 2, nrow = length(ct_mk_l_filt))
colnames(n_de) = c("SS", "1_wpi_neg")
rownames(n_de) = names(ct_mk_l_filt)
for(n in names(ct_mk_l_filt)){
tab = table(ct_mk_l_filt[[n]]$group)
n_de[n, "SS"] = tab["SS"]
n_de[n, "1_wpi_neg"] = tab["1_wpi_neg"]
}
n_de[is.na(n_de)] = 0
n_de_cl = matrix(0, ncol = 2, nrow = length(cl_mk_l_filt))
colnames(n_de_cl) = c("SS", "1_wpi_neg")
rownames(n_de_cl) = names(cl_mk_l_filt)
for(n in names(cl_mk_l_filt)){
tab = table(cl_mk_l_filt[[n]]$group)
n_de_cl[n, "SS"] = tab["SS"]
n_de_cl[n, "1_wpi_neg"] = tab["1_wpi_neg"]
}
n_de_cl[is.na(n_de_cl)] = 0
Compare clusters and conditions and cell types
df_ctcl = data.frame(table(ax_neg_srat$ctlabels, ax_neg_srat$integNeg_res.10))
pheatmap::pheatmap(table(ax_neg_srat$ctlabels, ax_neg_srat$integNeg_res.10),
fontsize = 6.5, clustering_method = "ward.D", scale = "row")
df_clct = data.frame(table(ax_neg_srat$integNeg_res.10, ax_neg_srat$sample_simp))
tab_s = prop.table(table(ax_neg_srat$integNeg_res.10, ax_neg_srat$sample_simp), margin = 2)
pheatmap::pheatmap(tab_s, fontsize = 6.5, clustering_method = "ward.D")
pheatmap::pheatmap(tab_s[tab_s[,1]>tab_s[,2]*1.1,], fontsize = 6.5, clustering_method = "ward.D")
Count number of DE genes - plot
Save data
plot_df = reshape2::melt(n_de)
ord = rowSums(n_de)
plot_df$Var1 = factor(plot_df$Var1, levels = rownames(n_de)[order(ord, decreasing = T)])
ggplot(plot_df, aes(x = Var1, y = value, fill = Var2))+
geom_col()+
scale_y_continuous(expand = c(0,0), limits = c(0, max(ord)+max(ord)*0.05))+
labs(x = "Cell Types", y = "# DE genes")+
theme_bw()+
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))